Evaluation of the impacts of regeneration on soil carbon and biodiversity in the Scottish Highlands
2025-03-2
Install Packages
#library(nlme)
#library(ggplot2)
#library(dplyr)
#library(broom)
#library(purrr)
#library(tidyr)
#library(readxl)
#library(RColorBrewer)
#library(gridExtra)
#library(vegan)
#library(patchwork)
#library(cowplot)
#library(grid)
#library(lme4)
#library(plotly)
#library(MuMIn)
#library(lmerTest)
#library(arm)
#library(sjPlot)
#library(indicspecies)
#library(stringr)Load Packages
library(nlme)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(broom)
library(purrr)
library(tidyr)
library(readxl)
library(RColorBrewer)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(vegan)## Loading required package: permute## Loading required package: latticelibrary(patchwork)
library(cowplot)##
## Attaching package: 'cowplot'## The following object is masked from 'package:patchwork':
##
## align_plotslibrary(grid)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(MuMIn)## Registered S3 methods overwritten by 'MuMIn':
## method from
## nobs.multinom broom
## nobs.fitdistr broomlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(arm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:plotly':
##
## select## The following object is masked from 'package:patchwork':
##
## area## The following object is masked from 'package:dplyr':
##
## select##
## arm (Version 1.14-4, built: 2024-4-1)## Working directory is /cloud/project##
## Attaching package: 'arm'## The following object is masked from 'package:MuMIn':
##
## coefplotlibrary(sjPlot)## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!##
## Attaching package: 'sjPlot'## The following objects are masked from 'package:cowplot':
##
## plot_grid, save_plotlibrary(indicspecies)
library(stringr)Functions
save_plot <- function(plot_obj, file_name, width = 8, height = 6, dpi = 300, format = "png", plot_title = NULL) {
format <- tolower(format)
# Update plot title if provided
if (!is.null(plot_title)) {
plot_obj <- plot_obj + ggtitle(plot_title)
}
# Choose graphics device
if (format == "png") {
png(paste0(file_name, ".png"), width = width, height = height, units = "in", res = dpi)
} else if (format == "pdf") {
pdf(paste0(file_name, ".pdf"), width = width, height = height)
} else if (format == "jpeg") {
jpeg(paste0(file_name, ".jpeg"), width = width, height = height, units = "in", res = dpi)
} else {
stop("Unsupported file format")
}
# Print and close
print(plot_obj)
dev.off()
}
make_seedling_plot <- function(data, group_label, color_line = "steelblue", color_point = "steelblue4", custom_title = NULL) {
# Fit model
model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha,
random = ~ 1 | Site,
data = data)
# Summary stats
model_summary <- summary(model)
p_val <- model_summary$tTable["Carbon_tonne_ha", "p-value"]
slope <- model_summary$tTable["Carbon_tonne_ha", "Value"]
r2_vals <- r.squaredGLMM(model)
r2_marginal <- r2_vals[1]
r2_conditional <- r2_vals[2]
label_text <- paste0(
"beta = ", round(slope, 3), "\n",
"R2[m] = ", round(r2_marginal, 3), "\n",
"R2[c] = ", round(r2_conditional, 3), "\n",
"p = ", signif(p_val, 3)
)
# Plot
graph <- plot_model(model, type = "pred", colors = color_line) +
geom_point(data = data, aes(x = Carbon_tonne_ha, y = sum_seedlings),
color = color_point, size = 1.5) +
labs(
x = "Carbon (tonne/ha)",
y = "Total Seedlings",
title = ifelse(is.null(custom_title), paste("Seedling Abundance vs topsoil SOC (", group_label, ")", sep = ""), custom_title)
) +
annotate("text", x = Inf, y = Inf, label = label_text,
hjust = 1.1, vjust = 1.5, size = 4, family = "mono") +
theme_bw() +
theme(
)
return(graph)
}Import Data
carbon_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 10)
overview_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 2)
seedling_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 6)
depth_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 9)## New names:
## • `` -> `...10`plant_com_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 3)
pitfall_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 8)
soil_invert_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 11)Process Data
#Mutate overview dataset to convert "Unforested (heath)", "Unforested (heathland)", "Unforested (heath)" and "Unforested (remnant_scat)" to "Unforested (remnant)"
overview_dataset <- overview_dataset %>%
mutate(Type = ifelse(Type %in% c("Unforested (heath)", "Unforested (heathland)"), "Unforested (heath)", Type)) %>%
mutate(Type = ifelse(Type %in% c("Unforested (remnant_scat)"), "Unforested (remnant)", Type))
carbon_dataset <- merge(carbon_dataset, overview_dataset, by = "Plot")
plant_com_dataset <- merge(plant_com_dataset, overview_dataset, by = "Plot")
plant_com_dataset <- plant_com_dataset %>%
mutate(Species = str_replace(Species, "Clubmoss", "Club moss"))
# Subset the data set to include only relevant information:
carbon_dataset <- carbon_dataset[, c("Site", "Plot","Layer" ,"Type", "Carbon (tonne/ha)")]
carbon_dataset <- carbon_dataset %>%
rename(Carbon_tonne_ha = `Carbon (tonne/ha)`)
# Subset the data set to include only relevant information:
carbon_dataset <- carbon_dataset[, c("Site", "Plot","Layer" ,"Type", "Carbon_tonne_ha")]
sum_seedlings <- seedling_dataset %>%
group_by(Plot) %>%
summarise(sum_seedlings = sum(Frequency), na.rm = TRUE)
est_seedlings <- seedling_dataset %>%
filter(Height_category %in% c("100-150", ">150", "50-100"))
est_seedlings <- est_seedlings %>%
group_by(Plot) %>% # Group by Plot
summarise(Total_Seedlings = sum(Frequency, na.rm = TRUE))
est_seedlings$sum_seedlings <- NULL # Remove the column
names(est_seedlings)[names(est_seedlings) == "Total_Seedlings"] <- "sum_seedlings" # Rename
carbon_dataset <- carbon_dataset %>%
filter(`Carbon_tonne_ha` != 0)
# First, get the plots with 0-10, O
o_samples <- carbon_dataset %>%
filter(Layer == "0-10, O")
# Then, get the plots with 0-10, M, but only those that do not have 0-10, O
m_samples <- carbon_dataset %>%
filter(Layer == "0-10, M") %>%
anti_join(o_samples, by = "Plot") # Exclude plots that have 0-10, O
# Combine both datasets (0-10, O + the filtered 0-10, M)
upper10_data <- bind_rows(o_samples, m_samples)
upper10_seedlings <- merge(upper10_data, sum_seedlings, by = "Plot")
upper10_est_seedlings <- merge(upper10_data, est_seedlings, by = "Plot")
regen_upper10 <- upper10_seedlings[upper10_seedlings$Type == "Regen", ]
regen_upper10_est <- upper10_est_seedlings[upper10_est_seedlings$Type == "Regen", ]
carbon_0_20 <- carbon_dataset %>%
# Create a new layer that combines "0-10, O" and "10-20, O"
mutate(Layer_combined = case_when(
Layer == "0-10, O" ~ "0-20, O",
Layer == "10-20, O" ~ "0-20, O",
TRUE ~ Layer # Keep other layers unchanged
)) %>%
# Now group by the combined layer and other factors (Type, Site)
group_by(Plot, Type, Layer_combined, Site) %>%
summarise(
Total_Carbon_Stock = sum(Carbon_tonne_ha) # Upper bound of 95% CI
)## `summarise()` has grouped output by 'Plot', 'Type', 'Layer_combined'. You can
## override using the `.groups` argument.plant_com_dataset <- plant_com_dataset %>%
mutate(DOMIN = case_when(
DOMIN == 1 ~ 1,
DOMIN == 2 ~ 2,
DOMIN == 3 ~ 3,
DOMIN == 4 ~ 7,
DOMIN == 5 ~ 18,
DOMIN == 6 ~ 29.5,
DOMIN == 7 ~ 42,
DOMIN == 8 ~ 63,
DOMIN == 9 ~ 83,
DOMIN == 10 ~ 95.5 # For any values outside 1-10
))
plant_com_dataset <- plant_com_dataset[c("Plot", "Site.y", "Species", "DOMIN", "Type")]
plant_rich <- plant_com_dataset %>%
group_by(Plot) %>%
summarise(Richness = sum(DOMIN > 0))
plant_rich <- merge(plant_rich, overview_dataset, by = "Plot")
# Calculate Simpson's Diversity Index for each plot
plant_div <- plant_com_dataset %>%
group_by(Plot) %>%
# Calculate total abundance for each plot
mutate(Total_DOMIN = sum(DOMIN)) %>%
# Calculate relative abundance for each species
mutate(Relative_DOMIN = DOMIN / Total_DOMIN) %>%
# Calculate Simpson's D (sum of squared relative abundances)
summarise(Simpsons_D = sum(Relative_DOMIN^2, na.rm = TRUE)) %>%
# Calculate Simpson's Diversity Index (1 - D)
mutate(Diversity_Index = 1 - Simpsons_D)
plant_div <- merge(plant_div, overview_dataset, by = "Plot")
soil_invert_dataset <- merge(soil_invert_dataset, overview_dataset, by = "Plot")
soil_invert_dataset <- soil_invert_dataset[c("Plot", "Site.y", "Order", "Abundance", "Type")]
soil_invert_rich <- soil_invert_dataset %>%
group_by(Plot) %>%
summarise(Richness = sum(Abundance > 0))
soil_invert_rich <- merge(soil_invert_rich, overview_dataset, by = "Plot")
# Calculate Simpson's Diversity Index for each plot
soil_invert_div <- soil_invert_dataset %>%
group_by(Plot) %>%
# Calculate total abundance for each plot
mutate(Total_Abundance = sum(Abundance)) %>%
# Calculate relative abundance for each species
mutate(Relative_Abundance = Abundance / Total_Abundance) %>%
# Calculate Simpson's D (sum of squared relative abundances)
summarise(Simpsons_D = sum(Relative_Abundance^2, na.rm = TRUE)) %>%
# Calculate Simpson's Diversity Index (1 - D)
mutate(Diversity_Index = 1 - Simpsons_D)
soil_invert_div <- merge(soil_invert_div, overview_dataset, by = "Plot")
pitfall_dataset <- merge(pitfall_dataset, overview_dataset, by = "Plot")
pitfall_dataset <- pitfall_dataset[c("Plot", "Site.y", "Order", "Abundance", "Type")]
pitfall_rich <- pitfall_dataset %>%
group_by(Plot) %>%
summarise(Richness = sum(Abundance > 0))
pitfall_rich <- pitfall_rich %>%
filter(Richness > 0)
pitfall_rich <- merge(pitfall_rich, overview_dataset, by = "Plot")
# Calculate Simpson's Diversity Index for each plot
pitfall_div <- pitfall_dataset %>%
group_by(Plot) %>%
# Calculate total abundance for each plot
mutate(Total_Abundance = sum(Abundance)) %>%
# Calculate relative abundance for each species
mutate(Relative_Abundance = Abundance / Total_Abundance) %>%
# Calculate Simpson's D (sum of squared relative abundances)
summarise(Simpsons_D = sum(Relative_Abundance^2, na.rm = TRUE)) %>%
# Calculate Simpson's Diversity Index (1 - D)
mutate(Diversity_Index = 1 - Simpsons_D)
pitfall_div <- merge(pitfall_div, overview_dataset, by = "Plot")
pitfall_div <- pitfall_div[pitfall_div$Diversity_Index != 0, ]Model Evaluation:
Soil Carbon/Woodland Type Model - Non-Normal residuals
### Soil Carbon/Woodland Type Model - Normal residuals
###############################################################################################
# Soil Carbon/Woodland Type Model:
soil_carbon_model <- lme(fixed = Total_Carbon_Stock ~ Type * Layer_combined,
random = ~ 1 | Site, # Random intercepts and slopes for Type within Site
data = carbon_0_20,
control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))
# Soil Carbon Model R-squared:
soil_carbon_r_squared <- r.squaredGLMM(soil_carbon_model)
# Post-Hoc Results:
summary(soil_carbon_model)## Linear mixed-effects model fit by REML
## Data: carbon_0_20
## AIC BIC logLik
## 1278.591 1306.713 -629.2954
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 14.93459 36.21726
##
## Fixed effects: Total_Carbon_Stock ~ Type * Layer_combined
## Value Std.Error DF
## (Intercept) 49.11022 12.54110 121
## TypeRegen -13.15502 12.61891 121
## TypeUnforested (heath) -9.50266 13.54644 121
## TypeUnforested (remnant) -11.34849 12.80473 121
## Layer_combined0-20, O 8.34979 13.09081 121
## TypeRegen:Layer_combined0-20, O 11.39713 17.68497 121
## TypeUnforested (heath):Layer_combined0-20, O -17.08628 19.10927 121
## TypeUnforested (remnant):Layer_combined0-20, O 3.81800 17.80741 121
## t-value p-value
## (Intercept) 3.915942 0.0001
## TypeRegen -1.042485 0.2993
## TypeUnforested (heath) -0.701488 0.4843
## TypeUnforested (remnant) -0.886273 0.3772
## Layer_combined0-20, O 0.637836 0.5248
## TypeRegen:Layer_combined0-20, O 0.644453 0.5205
## TypeUnforested (heath):Layer_combined0-20, O -0.894136 0.3730
## TypeUnforested (remnant):Layer_combined0-20, O 0.214405 0.8306
## Correlation:
## (Intr) TypRgn TypU(h) TypU(r)
## TypeRegen -0.520
## TypeUnforested (heath) -0.487 0.481
## TypeUnforested (remnant) -0.511 0.507 0.473
## Layer_combined0-20, O -0.507 0.499 0.469 0.489
## TypeRegen:Layer_combined0-20, O 0.372 -0.714 -0.344 -0.362
## TypeUnforested (heath):Layer_combined0-20, O 0.346 -0.341 -0.710 -0.335
## TypeUnforested (remnant):Layer_combined0-20, O 0.368 -0.365 -0.341 -0.719
## L_0-2O TR:L_O TypU(h):L_0-20,O
## TypeRegen
## TypeUnforested (heath)
## TypeUnforested (remnant)
## Layer_combined0-20, O
## TypeRegen:Layer_combined0-20, O -0.736
## TypeUnforested (heath):Layer_combined0-20, O -0.683 0.503
## TypeUnforested (remnant):Layer_combined0-20, O -0.729 0.539 0.499
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.2472249 -0.6413349 -0.1636584 0.4138244 5.5070856
##
## Number of Observations: 131
## Number of Groups: 3print(soil_carbon_r_squared)## R2m R2c
## [1,] 0.05122386 0.1891091###############################################################################################
# Soil Carbon/Woodland Type Model Residuals Normal Fit
residuals <- residuals(soil_carbon_model)
# Histogram to visualize residuals
hist(residuals, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black", breaks = 20)
# Q-Q plot to check for normality
qqnorm(residuals)
qqline(residuals)
# Shapiro-Wilk normality test
shapiro.test(residuals)##
## Shapiro-Wilk normality test
##
## data: residuals
## W = 0.80054, p-value = 4.619e-12###############################################################################################Seedling/Soil Carbon Model - Normal residuals
###############################################################################################
# Seedling/Soil Carbon Model:
upper10_seedling_model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha,
random = ~ 1 | Site,
data = regen_upper10)
upper10_est_seedling_model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha,
random = ~ 1 | Site,
data = regen_upper10_est)
# Seedling/Soil Carbon R-squared:
upper10_seedling_r_squared <- r.squaredGLMM(upper10_seedling_model)
upper10_est_seedling_r_squared <- r.squaredGLMM(upper10_est_seedling_model)
# Post-Hoc Results
summary(upper10_seedling_model)## Linear mixed-effects model fit by REML
## Data: regen_upper10
## AIC BIC logLik
## 209.4996 213.2773 -100.7498
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 23.32199 32.94415
##
## Fixed effects: sum_seedlings ~ Carbon_tonne_ha
## Value Std.Error DF t-value p-value
## (Intercept) 90.73903 20.533222 17 4.419132 0.0004
## Carbon_tonne_ha -0.51144 0.420827 17 -1.215323 0.2409
## Correlation:
## (Intr)
## Carbon_tonne_ha -0.669
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.99297649 -0.52684978 0.04632305 0.63364651 1.38446672
##
## Number of Observations: 21
## Number of Groups: 3print(upper10_seedling_r_squared)## R2m R2c
## [1,] 0.0578357 0.372375summary(upper10_est_seedling_model)## Linear mixed-effects model fit by REML
## Data: regen_upper10_est
## AIC BIC logLik
## 172.6892 176.022 -82.34459
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 26.43445 18.74382
##
## Fixed effects: sum_seedlings ~ Carbon_tonne_ha
## Value Std.Error DF t-value p-value
## (Intercept) 41.4879 17.738754 15 2.3388286 0.0336
## Carbon_tonne_ha -0.0868 0.246273 15 -0.3524518 0.7294
## Correlation:
## (Intr)
## Carbon_tonne_ha -0.447
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.69250319 -0.43396292 -0.03079282 0.60978544 1.48404159
##
## Number of Observations: 19
## Number of Groups: 3print(upper10_est_seedling_r_squared)## R2m R2c
## [1,] 0.002990246 0.666435###############################################################################################
# Seedling/Soil Carbon Model Residuals Normal Fit
residuals_upper10_seedlings <- residuals(upper10_seedling_model)
residuals_upper10_est_seedlings <- residuals(upper10_est_seedling_model)
# Histogram to visualize residuals
hist(residuals_upper10_seedlings, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Histogram to visualize residuals
hist(residuals_upper10_est_seedlings, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_upper10_seedlings)
qqline(residuals_upper10_seedlings)
# Shapiro-Wilk normality test
shapiro.test(residuals_upper10_seedlings)##
## Shapiro-Wilk normality test
##
## data: residuals_upper10_seedlings
## W = 0.96306, p-value = 0.5797# Q-Q plot to check for normality
qqnorm(residuals_upper10_est_seedlings)
qqline(residuals_upper10_est_seedlings)
# Shapiro-Wilk normality test
shapiro.test(residuals_upper10_est_seedlings)##
## Shapiro-Wilk normality test
##
## data: residuals_upper10_est_seedlings
## W = 0.91636, p-value = 0.09694###############################################################################################Plant Diversity Model - Non-Normal residuals
###############################################################################################
# Plant Diversity Model
diversity_model <- lme(fixed = Diversity_Index ~ Type,
random = ~ 1 | Site,
data = plant_div,
control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))
# Plant Diversity Model R-squared:
plant_diversity_r_squared <- r.squaredGLMM(diversity_model)
summary(diversity_model)## Linear mixed-effects model fit by REML
## Data: plant_div
## AIC BIC logLik
## -79.72206 -65.89767 45.86103
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 0.003486986 0.1201534
##
## Fixed effects: Diversity_Index ~ Type
## Value Std.Error DF t-value p-value
## (Intercept) 0.6861663 0.02629680 72 26.093143 0.0000
## TypeRegen -0.0356075 0.03708015 72 -0.960284 0.3401
## TypeUnforested (heath) 0.0490610 0.04061927 72 1.207826 0.2311
## TypeUnforested (remnant) 0.0139425 0.03708015 72 0.376009 0.7080
## Correlation:
## (Intr) TypRgn TypU(h)
## TypeRegen -0.705
## TypeUnforested (heath) -0.644 0.456
## TypeUnforested (remnant) -0.705 0.500 0.456
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.8956063 -0.5227355 0.1062503 0.6021995 1.6875788
##
## Number of Observations: 78
## Number of Groups: 3print(plant_diversity_r_squared)## R2m R2c
## [1,] 0.05584031 0.05663483###############################################################################################
# Plant Diversity Model Residuals Normal Fit
residuals_plant_diversity <- residuals(diversity_model)
# Histogram to visualize residuals
hist(residuals_plant_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_plant_diversity)
qqline(residuals_plant_diversity)
# Shapiro-Wilk normality test
shapiro.test(residuals_plant_diversity)##
## Shapiro-Wilk normality test
##
## data: residuals_plant_diversity
## W = 0.94101, p-value = 0.001287###############################################################################################Plant Richness Model - Not normal residuals
# Mixed model (if it fits)
richness_model_lme <- glmer(Richness ~ Type + (1 | Site), data = plant_rich, family = poisson)## boundary (singular) fit: see help('isSingular')# Fixed model
richness_model_lm <- lm(Richness ~ Type, data = plant_rich)
# Compare (will likely show negligible gain in model fit)
anova(richness_model_lme, richness_model_lm)## Data: plant_rich
## Models:
## richness_model_lme: Richness ~ Type + (1 | Site)
## richness_model_lm: Richness ~ Type
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## richness_model_lme 5 383.76 395.54 -186.88 373.76
## richness_model_lm 5 378.76 390.54 -184.38 368.76 5.002 0richness_model_lm <- glm(Richness ~ Type, data = plant_rich, family = poisson())
deviance <- deviance(richness_model_lm)
df_resid <- df.residual(richness_model_lm)
dispersion_ratio <- deviance / df_resid
dispersion_ratio## [1] 0.6040578###############################################################################################
# Plant Richness Model
richness_model <- glm(Richness ~ Type, data = plant_rich, family = poisson())
# Plant Diversity Model R-squared:
plant_richness_r_squared <- r.squaredGLMM(richness_model)## Warning: the null model is only correct if all the variables it uses are identical
## to those used in fitting the original model.summary(richness_model)##
## Call:
## glm(formula = Richness ~ Type, family = poisson(), data = plant_rich)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.339973 0.067729 34.549 <2e-16 ***
## TypeRegen 0.168465 0.091995 1.831 0.0671 .
## TypeUnforested (heath) 0.014572 0.104481 0.139 0.8891
## TypeUnforested (remnant) 0.009132 0.095565 0.096 0.9239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 49.301 on 77 degrees of freedom
## Residual deviance: 44.700 on 74 degrees of freedom
## AIC: 381.76
##
## Number of Fisher Scoring iterations: 4print(plant_richness_r_squared)## R2m R2c
## delta 0.05395323 0.05395323
## lognormal 0.05624468 0.05624468
## trigamma 0.05165377 0.05165377###############################################################################################
# Plant Diversity Model Residuals Normal Fit
residuals_plant_richness <- residuals(richness_model)
# Histogram to visualize residuals
hist(residuals_plant_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_plant_richness)
qqline(residuals_plant_richness)
# Shapiro-Wilk normality test
shapiro.test(residuals_plant_richness)##
## Shapiro-Wilk normality test
##
## data: residuals_plant_richness
## W = 0.93667, p-value = 0.0007625###############################################################################################Soil Invertebrate Diversity Model - Normal Residuals
###############################################################################################
# Soil Invertebrate Diversity Model
soil_invert_diversity_model <- lme(fixed = Diversity_Index ~ Type,
random = ~ 1 | Site,
data = soil_invert_div,
,
control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))
# Soil Invertebrate Diversity Model R-squared:
soil_invert_diversity_r_squared <- r.squaredGLMM(soil_invert_diversity_model)
summary(soil_invert_diversity_model)## Linear mixed-effects model fit by REML
## Data: soil_invert_div
## AIC BIC logLik
## -3.784492 7.187356 7.892246
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 0.00152749 0.1826833
##
## Fixed effects: Diversity_Index ~ Type
## Value Std.Error DF t-value p-value
## (Intercept) 0.3659797 0.04883611 45 7.494038 0.0000
## TypeRegen 0.0949750 0.07036316 45 1.349784 0.1838
## TypeUnforested (heath) -0.0654500 0.07563807 45 -0.865305 0.3915
## TypeUnforested (remnant) -0.0519503 0.07036316 45 -0.738317 0.4642
## Correlation:
## (Intr) TypRgn TypU(h)
## TypeRegen -0.694
## TypeUnforested (heath) -0.645 0.448
## TypeUnforested (remnant) -0.694 0.481 0.448
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.5234455 -0.6473544 0.1723194 0.5872532 1.9679246
##
## Number of Observations: 50
## Number of Groups: 2print(soil_invert_diversity_r_squared)## R2m R2c
## [1,] 0.106536 0.1065985###############################################################################################
# Soil Invertebrate Diversity Model Residuals Normal Fit
residuals_soil_invert_diversity <- residuals(soil_invert_diversity_model)
# Histogram to visualize residuals
hist(residuals_soil_invert_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_soil_invert_diversity)
qqline(residuals_soil_invert_diversity)
# Shapiro-Wilk normality test
shapiro.test(residuals_soil_invert_diversity)##
## Shapiro-Wilk normality test
##
## data: residuals_soil_invert_diversity
## W = 0.98845, p-value = 0.9036###############################################################################################Soil Invertebrate Richness Model - Not normal residuals
# Mixed model (if it fits)
soil_richness_invert_model_lme <- glmer(Richness ~ Type + (1 | Site), data = soil_invert_rich, family = poisson)## boundary (singular) fit: see help('isSingular')soil_invert_richness_model_lm <- glm(Richness ~ Type,
data = soil_invert_rich,
family = poisson)
# Compare (will likely show negligible gain in model fit)
anova(soil_richness_invert_model_lme, soil_invert_richness_model_lm)## Data: soil_invert_rich
## Models:
## soil_invert_richness_model_lm: Richness ~ Type
## soil_richness_invert_model_lme: Richness ~ Type + (1 | Site)
## npar AIC BIC logLik -2*log(L) Chisq Df
## soil_invert_richness_model_lm 4 168.72 176.37 -80.36 160.72
## soil_richness_invert_model_lme 5 170.72 180.28 -80.36 160.72 0 1
## Pr(>Chisq)
## soil_invert_richness_model_lm
## soil_richness_invert_model_lme 1deviance <- deviance(soil_invert_richness_model_lm)
df_resid <- df.residual(soil_invert_richness_model_lm)
dispersion_ratio <- deviance / df_resid
dispersion_ratio## [1] 0.3068677###############################################################################################
# Soil Invertebrate Richness Model
soil_invert_richness_model <- glm(Richness ~ Type,
data = soil_invert_rich,
family = poisson)
# Soil Invertebrate Richness Model R-squared:
soil_invert_richness_r_squared <- r.squaredGLMM(soil_invert_richness_model)## Warning: the null model is only correct if all the variables it uses are identical
## to those used in fitting the original model.summary(soil_invert_richness_model)##
## Call:
## glm(formula = Richness ~ Type, family = poisson, data = soil_invert_rich)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3312 0.1374 9.692 <2e-16 ***
## TypeRegen -0.2853 0.2142 -1.332 0.1830
## TypeUnforested (heath) -0.3380 0.2364 -1.429 0.1529
## TypeUnforested (remnant) -0.3997 0.2217 -1.802 0.0715 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 18.163 on 49 degrees of freedom
## Residual deviance: 14.116 on 46 degrees of freedom
## AIC: 168.72
##
## Number of Fisher Scoring iterations: 4print(soil_invert_richness_r_squared)## R2m R2c
## delta 0.07155917 0.07155917
## lognormal 0.08198368 0.08198368
## trigamma 0.06107931 0.06107931###############################################################################################
# Soil Invertebrate Richness Model Residuals Normal Fit
residuals_soil_invert_richness <- residuals(soil_invert_richness_model)
# Histogram to visualize residuals
hist(residuals_soil_invert_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_soil_invert_richness)
qqline(residuals_soil_invert_richness)
# Shapiro-Wilk normality test
shapiro.test(residuals_soil_invert_richness)##
## Shapiro-Wilk normality test
##
## data: residuals_soil_invert_richness
## W = 0.94179, p-value = 0.01587###############################################################################################Pitfall Invertebrate Diversity Model - Not Normal residuals
###############################################################################################
# Pitfall Invertebrate Diversity Model
pitfall_diversity_model <- lme(fixed = Diversity_Index ~ Type,
random = ~ 1 | Site,
data = pitfall_div,
,
control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))
# Soil Invertebrate Diversity Model R-squared:
pitfall_invert_diversity_r_squared <- r.squaredGLMM(pitfall_diversity_model)
summary(pitfall_diversity_model)## Linear mixed-effects model fit by REML
## Data: pitfall_div
## AIC BIC logLik
## -6.73401 6.756961 9.367005
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 0.05952625 0.1914478
##
## Fixed effects: Diversity_Index ~ Type
## Value Std.Error DF t-value p-value
## (Intercept) 0.6610537 0.05491415 68 12.037949 0.0000
## TypeRegen -0.0025667 0.06135071 68 -0.041836 0.9668
## TypeUnforested (heath) 0.0776850 0.06671556 68 1.164421 0.2483
## TypeUnforested (remnant) -0.0045245 0.05983130 68 -0.075620 0.9399
## Correlation:
## (Intr) TypRgn TypU(h)
## TypeRegen -0.543
## TypeUnforested (heath) -0.500 0.448
## TypeUnforested (remnant) -0.558 0.499 0.459
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.8213320 -0.4327584 0.2168915 0.6971038 1.6957427
##
## Number of Observations: 74
## Number of Groups: 3print(pitfall_invert_diversity_r_squared)## R2m R2c
## [1,] 0.02427668 0.1102898###############################################################################################
# Soil Invertebrate Diversity Model Residuals Normal Fit
residuals_pitfall_invert_diversity <- residuals(pitfall_diversity_model)
# Histogram to visualize residuals
hist(residuals_pitfall_invert_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_pitfall_invert_diversity)
qqline(residuals_pitfall_invert_diversity)
# Shapiro-Wilk normality test
shapiro.test(residuals_pitfall_invert_diversity)##
## Shapiro-Wilk normality test
##
## data: residuals_pitfall_invert_diversity
## W = 0.93568, p-value = 0.0009668###############################################################################################Pitfall Invertebrate Richness Model - Normal Residuals
###############################################################################################
# Pitfall Invertebrate Richness Model
pitfall_richness_model <- glmer(Richness ~ Type +
(1 | Site),
family = poisson(),
data = pitfall_rich)
# Soil Invertebrate Richness Model R-squared:
pitfall_invert_richness_r_squared <- r.squaredGLMM(pitfall_richness_model)## Warning: the null model is only correct if all the variables it uses are identical
## to those used in fitting the original model.summary(pitfall_richness_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Richness ~ Type + (1 | Site)
## Data: pitfall_rich
##
## AIC BIC logLik -2*log(L) df.resid
## 340.0 351.5 -165.0 330.0 69
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78501 -0.63764 -0.02627 0.67942 2.50383
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.02481 0.1575
## Number of obs: 74, groups: Site, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7718 0.1278 13.861 <2e-16 ***
## TypeRegen -0.2433 0.1348 -1.806 0.071 .
## TypeUnforested (heath) -0.1389 0.1540 -0.902 0.367
## TypeUnforested (remnant) -0.1735 0.1337 -1.298 0.194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TypRgn TypU(h)
## TypeRegen -0.463
## TypUnf(hth) -0.397 0.385
## TypU(rmnnt) -0.465 0.443 0.391print(pitfall_invert_richness_r_squared)## R2m R2c
## delta 0.03992091 0.1497831
## lognormal 0.04304071 0.1614886
## trigamma 0.03673198 0.1378183###############################################################################################
# Soil Invertebrate Richness Model Residuals Normal Fit
residuals_pitfall_invert_richness <- residuals(pitfall_richness_model)
# Histogram to visualize residuals
hist(residuals_pitfall_invert_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")
# Q-Q plot to check for normality
qqnorm(residuals_pitfall_invert_richness)
qqline(residuals_pitfall_invert_richness)
# Shapiro-Wilk normality test
shapiro.test(residuals_pitfall_invert_richness)##
## Shapiro-Wilk normality test
##
## data: residuals_pitfall_invert_richness
## W = 0.96675, p-value = 0.04836###############################################################################################Analysis:
Soil Carbon:
Extracting model predicted estimates:
library(nlme)
# Subset the data for Organic Soil (Layer_combined == '0-20, O')
organic_data <- carbon_0_20[carbon_0_20$Layer_combined == "0-20, O", ]
# Subset the data for Mineral Soil (Layer_combined == '0-10, M')
mineral_data <- carbon_0_20[carbon_0_20$Layer_combined == "0-10, M", ]
library(nlme)
organic_model <- lme(
Total_Carbon_Stock ~ Type,
random = ~ 1 | Site,
data = organic_data
)
mineral_model <- lme(
Total_Carbon_Stock ~ Type,
random = ~ 1 | Site,
data = mineral_data
)
organic_confints <- intervals(organic_model, level = 0.95)$fixed
mineral_confints <- intervals(mineral_model, which = "fixed", level = 0.95)
mineral_fixed <- mineral_confints$fixed
# Step 1: Get the intercept from the organic and mineral models
organic_intercept <- organic_confints["(Intercept)", "est."]
mineral_intercept <- mineral_fixed["(Intercept)", "est."] # or "Estimate" depending on version
# Step 2: Add the intercept to other treatment-level estimates
organic_adjusted <- organic_confints
organic_adjusted[-1, ] <- organic_adjusted[-1, ] + organic_intercept # Skip the intercept row
# Adjust the rest
mineral_adjusted <- mineral_fixed
mineral_adjusted[rownames(mineral_adjusted) != "(Intercept)", ] <-
mineral_adjusted[rownames(mineral_adjusted) != "(Intercept)", ] + mineral_intercept
mineral_adjusted## lower est. upper
## (Intercept) 36.97371 45.69073 54.40775
## TypeRegen 21.53700 33.68209 45.82717
## TypeUnforested (heath) 26.08765 39.10718 52.12671
## TypeUnforested (remnant) 22.01451 34.34224 46.66997
## attr(,"label")
## [1] "Fixed effects:"# Step 3: Rebuild the cleaned-up dataframes
organic_df <- data.frame(
Type = c("Mature", "Regen", "Unforested (heath)", "Unforested (remnant)"),
Soil_Type = "Organic",
Lower = organic_adjusted[, "lower"],
Estimate = organic_adjusted[, "est."],
Upper = organic_adjusted[, "upper"]
)
mineral_df <- data.frame(
Type = c("Mature", "Regen", "Unforested (heath)", "Unforested (remnant)"),
Soil_Type = "Mineral",
Lower = mineral_adjusted[, "lower"],
Estimate = mineral_adjusted[, "est."],
Upper = mineral_adjusted[, "upper"]
)
# Combine
combined_df <- rbind(organic_df, mineral_df)library(dplyr)
# Sort just to be safe
stacked_df <- combined_df %>%
arrange(Type, Soil_Type) %>%
mutate(
Soil_Type = factor(Soil_Type, levels = c("Organic", "Mineral"))
) %>%
group_by(Type) %>%
arrange(Soil_Type, .by_group = TRUE) %>%
mutate(
ymin = ifelse(Soil_Type == "Organic", 0, lag(Estimate)),
ymax = ymin + Estimate,
# For error bars
error_lower = ifelse(Soil_Type == "Organic", Lower, lag(Estimate) + Lower),
error_upper = ifelse(Soil_Type == "Organic", Upper, lag(Estimate) + Upper)
)
# Reorder the Soil_Type factor levels
stacked_df$Soil_Type <- factor(stacked_df$Soil_Type, levels = c("Mineral", "Organic")) # Change order here
stacked_df## # A tibble: 8 × 9
## # Groups: Type [4]
## Type Soil_Type Lower Estimate Upper ymin ymax error_lower error_upper
## <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Mature Organic 21.7 56.7 91.7 0 56.7 21.7 91.7
## 2 Mature Mineral 37.0 45.7 54.4 56.7 102. 93.7 111.
## 3 Regen Organic 24.0 55.6 87.1 0 55.6 24.0 87.1
## 4 Regen Mineral 21.5 33.7 45.8 55.6 89.2 77.1 101.
## 5 Unforested… Organic -3.72 30.6 64.9 0 30.6 -3.72 64.9
## 6 Unforested… Mineral 26.1 39.1 52.1 30.6 69.7 56.7 82.7
## 7 Unforested… Organic 18.2 49.7 81.2 0 49.7 18.2 81.2
## 8 Unforested… Mineral 22.0 34.3 46.7 49.7 84.1 71.7 96.4organic_estimates <- stacked_df %>%
dplyr::filter(Soil_Type == "Organic") %>% # Filter only organic estimates
dplyr::select(Type, Estimate) # Select Type and Estimate columns
# View the organic estimates
print(organic_estimates)## # A tibble: 4 × 2
## # Groups: Type [4]
## Type Estimate
## <chr> <dbl>
## 1 Mature 56.7
## 2 Regen 55.6
## 3 Unforested (heath) 30.6
## 4 Unforested (remnant) 49.7# Join the organic estimates with the raw carbon data
carbon_0_20_adjusted <- carbon_0_20 %>%
left_join(organic_estimates, by = "Type") # Join by 'Type' (treatment)
# View the adjusted data to confirm the join
head(carbon_0_20_adjusted)## # A tibble: 6 × 6
## # Groups: Plot, Type, Layer_combined [6]
## Plot Type Layer_combined Site Total_Carbon_Stock Estimate
## <dbl> <chr> <chr> <chr> <dbl> <dbl>
## 1 1 Unforested (remnant) 0-20, O Glen Ma… 24.5 49.7
## 2 2 Unforested (remnant) 0-20, O Glen Ma… 31.2 49.7
## 3 3 Mature 0-20, O Glen Ma… 193. 56.7
## 4 4 Mature 0-10, M Glen Ma… 85.0 56.7
## 5 5 Unforested (remnant) 0-10, M Glen Ma… 52.3 49.7
## 6 5 Unforested (remnant) 0-20, O Glen Ma… 72.3 49.7# Adjust the mineral carbon stock values by adding the organic estimate
carbon_0_20_adjusted <- carbon_0_20_adjusted %>%
mutate(
Total_Carbon_Stock_adjusted = ifelse(
Layer_combined == 'Mineral', # If it's mineral soil
Total_Carbon_Stock + Estimate, # Add the organic estimate
Total_Carbon_Stock # Keep the original carbon stock for organic soil
)
)
# View the adjusted data
head(carbon_0_20_adjusted)## # A tibble: 6 × 7
## # Groups: Plot, Type, Layer_combined [6]
## Plot Type Layer_combined Site Total_Carbon_Stock Estimate
## <dbl> <chr> <chr> <chr> <dbl> <dbl>
## 1 1 Unforested (remnant) 0-20, O Glen Ma… 24.5 49.7
## 2 2 Unforested (remnant) 0-20, O Glen Ma… 31.2 49.7
## 3 3 Mature 0-20, O Glen Ma… 193. 56.7
## 4 4 Mature 0-10, M Glen Ma… 85.0 56.7
## 5 5 Unforested (remnant) 0-10, M Glen Ma… 52.3 49.7
## 6 5 Unforested (remnant) 0-20, O Glen Ma… 72.3 49.7
## # ℹ 1 more variable: Total_Carbon_Stock_adjusted <dbl># Recode the 'Layer_combined' column to use 'Mineral' and 'Organic'
carbon_0_20_adjusted <- carbon_0_20_adjusted %>%
mutate(
Layer_combined = recode(Layer_combined,
'0-10, M' = 'Mineral', # Change '0-10, M' to 'Mineral'
'0-20, O' = 'Organic') # Optionally, recode 'Organic' if needed
)Plotting soil carbon stock
stacked_df$Treatment_Layer <- paste(stacked_df$Type, stacked_df$Soil_Type, sep = " / ")
carbon_0_20_adjusted$Treatment_Layer <- paste(carbon_0_20_adjusted$Type, carbon_0_20_adjusted$Layer_combined, sep = " / ")
library(colorspace)
palette.colors()## [1] "#000000" "#E69F00" "#56B4E9" "#009E73" "#F0E442" "#0072B2" "#D55E00"
## [8] "#CC79A7" "#999999"# Base woodland treatment colors
woodland_base <- c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)
# Build combined color palette
layer_colors <- c()
for (type in names(woodland_base)) {
layer_colors[paste0(type, " / Organic")] <- woodland_base[type]
layer_colors[paste0(type, " / Mineral")] <- darken(woodland_base[type], amount = 0.2)
}
library(ggplot2)
soil_woodland_graph <- ggplot(stacked_df, aes(x = Type, y = Estimate, fill = Treatment_Layer)) +
geom_bar(stat = "identity", position = "stack", color = "black", size = 0.5, width = 0.6, alpha = 0.7) +
geom_errorbar(aes(ymin = error_lower, ymax = error_upper),
position = "identity", width = 0.3, color = "black", size = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
scale_y_reverse() +
scale_shape_manual(values = c("Organic" = 21, "Mineral" = 24)) +
scale_fill_manual(values = layer_colors) +
guides(
fill = guide_legend(override.aes = list(shape = c(21, 24)), reverse = TRUE),
shape = "none" # Hide the separate shape legend
) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1), # <-- added line
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10)
) +
labs(
x = "Treatment Type",
y = "Carbon Stock (tonne/ha)",
fill = "Forest Treatment / Soil Layer",
title = "Stacked Organic & Mineral Soil Carbon by Treatment"
) +
geom_jitter(data = carbon_0_20_adjusted,
aes(x = Type,
y = ifelse(Layer_combined == 'Mineral',
Total_Carbon_Stock_adjusted,
Total_Carbon_Stock),
fill = Treatment_Layer,
shape = Layer_combined),
width = 0.15,
height = 0,
size = 1.5,
color = "black",
stroke = 0.8,
alpha = 0.5)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.soil_woodland_graph
save_plot(soil_woodland_graph, "soil_carbon_graph", plot_title = "Carbon stock by Forest Treatment")## png
## 2Plotting Organic Depth Horizon
# Modify the O_layer_depth column
depth_dataset <- depth_dataset %>%
mutate(
O_layer_depth = case_when(
O_layer_depth == ">20" ~ 25, # Replace ">20" with 25
O_layer_depth == "(=20)" ~ 20, # Replace "(=20)" with 20
TRUE ~ as.numeric(O_layer_depth) # Convert remaining values to numeric
)
) %>%
filter(!is.na(O_layer_depth)) # Remove rows where O_layer_depth is NA## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `O_layer_depth = case_when(...)`.
## Caused by warning:
## ! NAs introduced by coerciondepth_dataset <- depth_dataset[, c("Plot", "O_layer_depth", "Site")]
depth_dataset <- merge(depth_dataset, overview_dataset, by = c("Plot", "Site"))
depth_dataset <- depth_dataset[, c("Plot", "O_layer_depth", "Type", "Site")]
average_depth_results <- depth_dataset %>%
group_by(Type, Site) %>%
summarise(
Average_O_layer_depth = mean(O_layer_depth, na.rm = TRUE), # Mean depth
SE_O_layer_depth = sd(O_layer_depth, na.rm = TRUE) / sqrt(n()), # Standard Error
Lower_CI_depth = Average_O_layer_depth - 1.96 * SE_O_layer_depth, # Lower 95% CI
Upper_CI_depth = Average_O_layer_depth + 1.96 * SE_O_layer_depth # Upper 95% CI
)## `summarise()` has grouped output by 'Type'. You can override using the
## `.groups` argument.average_plot_depth <- depth_dataset %>%
group_by(Plot, Type, Site) %>%
summarise(
Plot_O_layer_depth = mean(O_layer_depth, na.rm = TRUE))## `summarise()` has grouped output by 'Plot', 'Type'. You can override using the
## `.groups` argument.woodland_base <- c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)
depth_graph <- ggplot(average_depth_results, aes(x = Type, y = Average_O_layer_depth, fill = Type)) +
geom_bar(stat = "identity", position = "stack", color = "black", size = 0.5, width = 0.6, alpha = 0.7) +
scale_fill_manual(name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
geom_errorbar(data = average_depth_results,
aes(ymin = Lower_CI_depth, ymax = Upper_CI_depth, x = Type),
width = 0.2, color = "black", size = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
scale_y_reverse() +
theme_minimal() + # Use a minimal theme with a white background
theme(
axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
strip.background = element_rect(fill = "grey", color = "black", size = 1), # Grey background for facet labels
strip.text = element_text(color = "black", size = 12), # Text color inside the grey box
panel.spacing = unit(1, "lines"), # Add spacing between panels if needed
panel.background = element_rect(fill = "white", color = "black"), # Apply white background with black border to panels
plot.margin = margin(10, 10, 10, 10) # Add margins around the plot if necessary
) +# Reverse y-axis for better visibility
labs(
x = "Forest Treatment", # Custom x-axis label
y = "Organic Layer Depth (cm)", # Custom y-axis label
fill = "Soil Layer", # Custom legend title
title = "Organic Horizon Depth across Sites" # Overall title
) + # Create facets based on 'Site'
theme(
axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_wrap(~Site)
depth_graph
save_plot(depth_graph, "depth_graph")## png
## 2# Perform one-way ANOVA for each site
anova_results <- depth_dataset %>%
do({
model <- aov(O_layer_depth ~ Type, data = .) # Perform ANOVA within each site
tidy(model) # Use the tidy() function from broom to get results in a nice format
})
# View the results of the ANOVA for each site
anova_results## # A tibble: 2 × 6
## term df sumsq meansq statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Type 3 429. 143. 1.58 0.195
## 2 Residuals 351 31836. 90.7 NA NA# Perform Tukey's HSD test for each site where ANOVA was significant
tukey_results <- depth_dataset %>%
group_by(Site) %>%
do({
model <- aov(O_layer_depth ~ Type, data = .) # Fit the ANOVA model
tukey <- TukeyHSD(model) # Perform Tukey's HSD post-hoc test
tidy(tukey) # Get results in a tidy format
})
# View Tukey HSD results for each site
tukey_results## # A tibble: 18 × 8
## # Groups: Site [3]
## Site term contrast null.value estimate conf.low conf.high adj.p.value
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Glen Cluan… Type Regen-M… 0 2.26 -3.01 7.53 0.680
## 2 Glen Cluan… Type Unfores… 0 1.22 -4.43 6.87 0.943
## 3 Glen Cluan… Type Unfores… 0 3.23 -2.00 8.46 0.377
## 4 Glen Cluan… Type Unfores… 0 -1.04 -6.53 4.45 0.960
## 5 Glen Cluan… Type Unfores… 0 0.974 -4.08 6.03 0.958
## 6 Glen Cluan… Type Unfores… 0 2.01 -3.44 7.47 0.771
## 7 Glen Loyne Type Regen-M… 0 1.26 -3.33 5.86 0.891
## 8 Glen Loyne Type Unfores… 0 1.93 -3.05 6.91 0.745
## 9 Glen Loyne Type Unfores… 0 1.25 -3.31 5.82 0.891
## 10 Glen Loyne Type Unfores… 0 0.666 -4.25 5.58 0.985
## 11 Glen Loyne Type Unfores… 0 -0.00840 -4.50 4.48 1.00
## 12 Glen Loyne Type Unfores… 0 -0.674 -5.56 4.21 0.984
## 13 Glen Mallie Type Regen-M… 0 0.866 -5.13 6.86 0.982
## 14 Glen Mallie Type Unfores… 0 -6.5 -13.7 0.657 0.0891
## 15 Glen Mallie Type Unfores… 0 0.883 -5.01 6.77 0.980
## 16 Glen Mallie Type Unfores… 0 -7.37 -14.6 -0.0871 0.0462
## 17 Glen Mallie Type Unfores… 0 0.0167 -6.02 6.06 1.00
## 18 Glen Mallie Type Unfores… 0 7.38 0.188 14.6 0.0420Seedling Establishmet on soil carbon:
Extracting and plotting model predicted data:
regen_upper10_graph <- make_seedling_plot(regen_upper10, "Regenerating", color_line = "#009E73", color_point = "black", custom_title = "Total Seedlings vs Topsoil SOC")
regen_upper10_est_graph <- make_seedling_plot(regen_upper10_est, "Regenerating", color_line = "#009E73", color_point = "black", custom_title = "Established Seedlings vs Topsoil SOC")
combined_regen_graph <- regen_upper10_graph / regen_upper10_est_graph
combined_regen_graph
save_plot(combined_regen_graph, "seedling_regen_graph")## png
## 2Plant Diversity:
Extracting model predicted data
predicted_plant_div_with_ci <- intervals(diversity_model, level = 0.95)
# Extract just the fixed effects into a data frame
plant_div_ci_df <- as.data.frame(predicted_plant_div_with_ci$fixed)
# Optionally, name the columns
colnames(plant_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")
# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
plant_div_ci_df$Effect <- rownames(predicted_plant_div_with_ci$fixed)
# Reorder for readability
plant_div_ci_df <- plant_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]
# View it
print(plant_div_ci_df)## Effect Estimate CI_lower
## (Intercept) (Intercept) 0.68616626 0.63374454
## TypeRegen TypeRegen -0.03560748 -0.10952542
## TypeUnforested (heath) TypeUnforested (heath) 0.04906101 -0.03191202
## TypeUnforested (remnant) TypeUnforested (remnant) 0.01394246 -0.05997548
## CI_upper
## (Intercept) 0.73858798
## TypeRegen 0.03831045
## TypeUnforested (heath) 0.13003405
## TypeUnforested (remnant) 0.08786039# Extract the intercept value
intercept_estimate <- plant_div_ci_df$Estimate[plant_div_ci_df$Effect == "(Intercept)"]
intercept_lower <- plant_div_ci_df$CI_lower[plant_div_ci_df$Effect == "(Intercept)"]
intercept_upper <- plant_div_ci_df$CI_upper[plant_div_ci_df$Effect == "(Intercept)"]
# Now create a new column with actual predicted values
plant_div_ci_df$Actual_Estimate <- plant_div_ci_df$Estimate
plant_div_ci_df$Actual_CI_lower <- plant_div_ci_df$CI_lower
plant_div_ci_df$Actual_CI_upper <- plant_div_ci_df$CI_upper
# Add the intercept to all *non-intercept* rows
non_intercept_rows <- plant_div_ci_df$Effect != "(Intercept)"
plant_div_ci_df$Actual_Estimate[non_intercept_rows] <-
plant_div_ci_df$Estimate[non_intercept_rows] + intercept_estimate
plant_div_ci_df$Actual_CI_lower[non_intercept_rows] <-
plant_div_ci_df$CI_lower[non_intercept_rows] + intercept_lower
plant_div_ci_df$Actual_CI_upper[non_intercept_rows] <-
plant_div_ci_df$CI_upper[non_intercept_rows] + intercept_upper
# Using base R
colnames(plant_div_ci_df)[colnames(plant_div_ci_df) == "Effect"] <- "Type"
plant_div_ci_df$Type[plant_div_ci_df$Type == "(Intercept)"] <- "Mature"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeRegen"] <- "Regen"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
summary(diversity_model)## Linear mixed-effects model fit by REML
## Data: plant_div
## AIC BIC logLik
## -79.72206 -65.89767 45.86103
##
## Random effects:
## Formula: ~1 | Site
## (Intercept) Residual
## StdDev: 0.003486986 0.1201534
##
## Fixed effects: Diversity_Index ~ Type
## Value Std.Error DF t-value p-value
## (Intercept) 0.6861663 0.02629680 72 26.093143 0.0000
## TypeRegen -0.0356075 0.03708015 72 -0.960284 0.3401
## TypeUnforested (heath) 0.0490610 0.04061927 72 1.207826 0.2311
## TypeUnforested (remnant) 0.0139425 0.03708015 72 0.376009 0.7080
## Correlation:
## (Intr) TypRgn TypU(h)
## TypeRegen -0.705
## TypeUnforested (heath) -0.644 0.456
## TypeUnforested (remnant) -0.705 0.500 0.456
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.8956063 -0.5227355 0.1062503 0.6021995 1.6875788
##
## Number of Observations: 78
## Number of Groups: 3plant_div_ci_df## Type Estimate CI_lower
## (Intercept) Mature 0.68616626 0.63374454
## TypeRegen Regen -0.03560748 -0.10952542
## TypeUnforested (heath) Unforested (heath) 0.04906101 -0.03191202
## TypeUnforested (remnant) Unforested (remnant) 0.01394246 -0.05997548
## CI_upper Actual_Estimate Actual_CI_lower
## (Intercept) 0.73858798 0.6861663 0.6337445
## TypeRegen 0.03831045 0.6505588 0.5242191
## TypeUnforested (heath) 0.13003405 0.7352273 0.6018325
## TypeUnforested (remnant) 0.08786039 0.7001087 0.5737691
## Actual_CI_upper
## (Intercept) 0.7385880
## TypeRegen 0.7768984
## TypeUnforested (heath) 0.8686220
## TypeUnforested (remnant) 0.8264484Plotting plant iversity by forest treatment
plant_diversity_graph <- ggplot(plant_div, aes(x = Type, y = Diversity_Index, fill = Type)) +
# Violin plot to show distribution of raw values
geom_violin(alpha = 0.3, color = NA, width = 0.9) +
# Individual data points (jittered)
geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
# Model estimates (mean points from plant_div_ci_df)
geom_point(data = plant_div_ci_df,
aes(x = Type, y = Actual_Estimate),
shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +
# Confidence intervals (error bars)
geom_errorbar(data = plant_div_ci_df,
aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper),
width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +
# Aesthetics and themes
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10),
legend.position = "none"
) +
labs(
x = "Forest Treatment",
y = "Diversity Index (1 - Simpson's D)",
fill = "Type",
title = "Plant Diversity"
) +
scale_fill_manual(values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))Plant Richness:
Extracting model predicted data:
# Step 1: Extract the intercept value
intercept_value <- coef(summary(richness_model))["(Intercept)", "Estimate"]
# Step 2: Get confidence intervals (for all fixed effects)
wald_ci <- confint(richness_model, level = 0.95)## Waiting for profiling to be done...# Step 3: Create results dataframe
results_df <- data.frame(
Effect = rownames(coef(summary(richness_model))),
Estimate = coef(summary(richness_model))[, "Estimate"],
Std_Error = coef(summary(richness_model))[, "Std. Error"],
p_value = coef(summary(richness_model))[, "Pr(>|z|)"],
CI_Lower = wald_ci[, 1],
CI_Upper = wald_ci[, 2]
)
# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
results_df$Real_Estimate <- results_df$Estimate
results_df$Real_CI_Lower <- results_df$CI_Lower
results_df$Real_CI_Upper <- results_df$CI_Upper
non_intercepts <- results_df$Effect != "(Intercept)"
results_df$Real_Estimate[non_intercepts] <- results_df$Estimate[non_intercepts] + intercept_value
results_df$Real_CI_Lower[non_intercepts] <- results_df$CI_Lower[non_intercepts] + intercept_value
results_df$Real_CI_Upper[non_intercepts] <- results_df$CI_Upper[non_intercepts] + intercept_value
# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(results_df)[colnames(results_df) == "Effect"] <- "Type"
results_df$Type[results_df$Type == "(Intercept)"] <- "Mature"
results_df$Type[results_df$Type == "TypeRegen"] <- "Regen"
results_df$Type[results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
results_df$Type[results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(results_df) <- results_df$Type
# Step 8: Exponentiate Estimates and Confidence Intervals
results_df$Real_Estimate <- exp(results_df$Real_Estimate)
results_df$Real_CI_Lower <- exp(results_df$Real_CI_Lower)
results_df$Real_CI_Upper <- exp(results_df$Real_CI_Upper)
print(results_df)## Type Estimate Std_Error p_value
## Mature Mature 2.339972625 0.06772851 1.460262e-261
## Regen Regen 0.168464522 0.09199522 6.706602e-02
## Unforested (heath) Unforested (heath) 0.014572207 0.10448093 8.890769e-01
## Unforested (remnant) Unforested (remnant) 0.009132484 0.09556467 9.238673e-01
## CI_Lower CI_Upper Real_Estimate Real_CI_Lower
## Mature 2.20422005 2.4698410 10.38095 9.063180
## Regen -0.01150943 0.3493513 12.28571 10.262158
## Unforested (heath) -0.19150620 0.2184142 10.53333 8.571703
## Unforested (remnant) -0.17828169 0.1966002 10.47619 8.685812
## Real_CI_Upper
## Mature 11.82057
## Regen 14.72172
## Unforested (heath) 12.91497
## Unforested (remnant) 12.63629Plotting plant richness by forest treatment
plant_richness_graph <- ggplot(plant_rich, aes(x = Type, y = Richness, fill = Type)) +
# Violin plot to show data distribution
geom_violin(alpha = 0.3, color = NA, width = 0.9) +
# Jittered individual data points
geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
# Model estimates (points)
geom_point(data = results_df, aes(x = Type, y = Real_Estimate),
shape = 21, fill = "black", size = 3, inherit.aes = FALSE) +
# Confidence intervals from model
geom_errorbar(data = results_df,
aes(x = Type, ymin = Real_CI_Lower, ymax = Real_CI_Upper),
width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10)
) +
labs(
x = "Forest Treatment",
y = "Species Richness",
fill = "Type",
title = "Plant Richness"
) +
scale_fill_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))
plant_rich_div_graph <- plant_diversity_graph + plant_richness_graph
plant_rich_div_graph
save_plot(plant_rich_div_graph, "plant_rich_div_graph")## png
## 2Plant Community Composition:
plant_com_dataset <- plant_com_dataset[-c(561, 701), ]
plant_com_wide <- as.data.frame(spread(plant_com_dataset, Species, DOMIN))After converting the data into wide format, I need to replace all the NA values with 0
plant_com_wide[is.na(plant_com_wide)] <- 0
plant_com_env <- plant_com_wide[c("Plot", "Site.y", "Type")]
plant_com_wide$Plot <- NULL
plant_com_wide$Site.y <- NULL
plant_com_wide$Type <- NULL
plant_com_wide$Family <- NULL
plant_com_wide$Genus <- NULL
plant_com_wide <- as.data.frame(lapply(plant_com_wide, as.integer))
str(plant_com_wide)
NMDS_plant_com <- metaMDS(plant_com_wide)
NMDS_plant_com
NMDS_plant_com$stress
stressplot(NMDS_plant_com)
plant_com_env$Type <- factor(plant_com_env$Type)
scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_plant_com, type = "n")
with(plant_com_env, points(NMDS_plant_com, display = "sites", col = colvec[Type],
scaling = scl, pch = 21, bg = colvec[Type]))
with(plant_com_env, legend("topleft", legend = levels(Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec))
NMDS_plant_com.fit <- envfit(NMDS_plant_com ~ Type , data = plant_com_env, perm = 999)
NMDS_plant_com.fit##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## TypeMature 0.8232 0.1681
## TypeRegen 0.0537 -0.1053
## TypeUnforested (heath) -0.6063 -0.0011
## TypeUnforested (remnant) -0.4438 -0.0620
##
## Goodness of fit:
## r2 Pr(>r)
## Type 0.5383 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(NMDS_plant_com.fit, dis = "site")
plot(NMDS_plant_com.fit)
### Plant NMDS
# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)
ordi_ellipses <- ordiellipse(NMDS_plant_com, plant_com_env$Type, kind = "se", conf = 0.95, draw = "none")
# 2. Extract ellipse coordinates into a data frame
ellipse_df <- do.call(rbind, lapply(names(ordi_ellipses), function(group) {
ellipse_points <- vegan:::veganCovEllipse(ordi_ellipses[[group]]$cov,
ordi_ellipses[[group]]$center,
ordi_ellipses[[group]]$scale)
data.frame(ellipse_points, Type = group)
}))
colnames(ellipse_df)[1:2] <- c("NMDS1", "NMDS2")
library(ggplot2)
# Your main NMDS points
nmds_df <- data.frame(
NMDS1 = NMDS_plant_com$points[, 1],
NMDS2 = NMDS_plant_com$points[, 2],
Type = plant_com_env$Type
)
anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(plant_com_wide)), plant_com_env$Type, distance = "bray", permutations = 9999)
R_val <- round(anosim_plant_com_r_tr$statistic, 3)
p_val <- ifelse(anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(anosim_plant_com_r_tr$signif, 2)))
anosim_label <- paste0("ANOSIM R = ", R_val, "\n", "p ", p_val)
# Site scores
site_scores <- as.data.frame(NMDS_plant_com$points)
colnames(site_scores)[1:2] <- c("NMDS1", "NMDS2")
site_scores$Type <- plant_com_env$Type # Add grouping
# Species scores (usually in $species, or use scores())
species_scores <- as.data.frame(scores(NMDS_plant_com, display = "species"))
species_scores$Species <- rownames(species_scores)
plant_nmds_hull_graph <- ggplot(nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
geom_point(shape = 21, size = 1, stroke = 1) +
geom_polygon(data = ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
alpha = 0.4, linewidth = 0.5) +
# ⬇️ Add species labels with inherit.aes = FALSE
geom_text(data = species_scores,
mapping = aes(x = NMDS1, y = NMDS2, label = Species),
inherit.aes = FALSE,
color = "white", size = 2, alpha = 0) +
scale_color_manual(name = "Forest Treatment",
values = colvec) +
scale_fill_manual(name = "Forest Treatment",
values = colvec) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10),
legend.position = "none"
) +
annotate("text", x = Inf, y = Inf, label = anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
labs(title = "Plant Community Separation");
plant_nmds_hull_graph <- plant_nmds_hull_graph + coord_fixed()
plant_nmds_hull_graph
save_plot(plant_nmds_hull_graph, "plant_nmds_hull_graph")## png
## 2Plant NMDS resampling
# Base ANOSIM R value with all species
base_R <- anosim(wisconsin(sqrt(plant_com_wide)), plant_com_env$Type, distance = "bray", permutations = 999)$statistic
# Store results
jackknife_results <- data.frame(Species = character(),
R_value = numeric(),
Delta_R = numeric(),
stringsAsFactors = FALSE)
# Loop over all species (columns in your wide matrix)
species_names <- colnames(plant_com_wide)
for (sp in species_names) {
# Remove one species
reduced_matrix <- plant_com_wide[, setdiff(species_names, sp)]
# Run ANOSIM
anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), plant_com_env$Type, distance = "bray", permutations = 999)
# Store R and change from base
R_val <- anosim_result$statistic
delta_R <- base_R - R_val
jackknife_results <- rbind(jackknife_results, data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}
# Sort by greatest impact (largest Delta_R)
jackknife_results <- jackknife_results[order(-jackknife_results$Delta_R), ]
plant_jk_graph <- ggplot(jackknife_results, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
geom_col(fill = "#0072B2") +
coord_flip() +
labs(title = "Plant Species Impact on NMDS Clustering",
x = "Species Removed",
y = "Decrease in ANOSIM R") +
theme_minimal()
plant_jk_graph
save_plot(plant_jk_graph, "plant_jk")## png
## 2Plant NMD (No Trees)
# Define species to exclude
tree_species <- c("Silver birch", "Rowan", "Downy birch")
# Filter matrix to exclude these species
plant_com_no_trees <- plant_com_wide[, !colnames(plant_com_wide) %in% tree_species]# Run NMDS on the modified dataset
NMDS_plant_com_no_trees <- metaMDS(wisconsin(sqrt(plant_com_no_trees)), distance = "bray", k = 2, trymax = 100)## Run 0 stress 0.1901061
## Run 1 stress 0.21296
## Run 2 stress 0.1883802
## ... New best solution
## ... Procrustes: rmse 0.03239242 max resid 0.1414564
## Run 3 stress 0.1890093
## Run 4 stress 0.194464
## Run 5 stress 0.1930042
## Run 6 stress 0.2142864
## Run 7 stress 0.2104666
## Run 8 stress 0.214057
## Run 9 stress 0.2076797
## Run 10 stress 0.1977019
## Run 11 stress 0.2114374
## Run 12 stress 0.2116637
## Run 13 stress 0.2044696
## Run 14 stress 0.196647
## Run 15 stress 0.1983483
## Run 16 stress 0.188837
## ... Procrustes: rmse 0.01125132 max resid 0.08060584
## Run 17 stress 0.2083896
## Run 18 stress 0.2083942
## Run 19 stress 0.2040663
## Run 20 stress 0.2083267
## Run 21 stress 0.1881442
## ... New best solution
## ... Procrustes: rmse 0.02078665 max resid 0.1242391
## Run 22 stress 0.2101414
## Run 23 stress 0.2157395
## Run 24 stress 0.2065798
## Run 25 stress 0.2136877
## Run 26 stress 0.1972916
## Run 27 stress 0.197115
## Run 28 stress 0.1943633
## Run 29 stress 0.2093802
## Run 30 stress 0.1890592
## Run 31 stress 0.1986611
## Run 32 stress 0.2160472
## Run 33 stress 0.2086404
## Run 34 stress 0.2102971
## Run 35 stress 0.1887508
## Run 36 stress 0.2195215
## Run 37 stress 0.1903416
## Run 38 stress 0.188585
## ... Procrustes: rmse 0.01487024 max resid 0.1208851
## Run 39 stress 0.1971506
## Run 40 stress 0.2079329
## Run 41 stress 0.1928808
## Run 42 stress 0.2134388
## Run 43 stress 0.219096
## Run 44 stress 0.2132646
## Run 45 stress 0.2222949
## Run 46 stress 0.198661
## Run 47 stress 0.2145521
## Run 48 stress 0.215177
## Run 49 stress 0.214961
## Run 50 stress 0.1881211
## ... New best solution
## ... Procrustes: rmse 0.001937691 max resid 0.01359463
## Run 51 stress 0.1900061
## Run 52 stress 0.2081687
## Run 53 stress 0.2173113
## Run 54 stress 0.1974102
## Run 55 stress 0.2001864
## Run 56 stress 0.2044841
## Run 57 stress 0.1961687
## Run 58 stress 0.2200975
## Run 59 stress 0.2147783
## Run 60 stress 0.2109233
## Run 61 stress 0.2042485
## Run 62 stress 0.2134556
## Run 63 stress 0.1941284
## Run 64 stress 0.1925906
## Run 65 stress 0.2060116
## Run 66 stress 0.1886733
## Run 67 stress 0.189057
## Run 68 stress 0.2013378
## Run 69 stress 0.2017538
## Run 70 stress 0.1881142
## ... New best solution
## ... Procrustes: rmse 0.002560617 max resid 0.01920495
## Run 71 stress 0.1885752
## ... Procrustes: rmse 0.01453685 max resid 0.1212143
## Run 72 stress 0.2116655
## Run 73 stress 0.2082069
## Run 74 stress 0.2039048
## Run 75 stress 0.1888904
## Run 76 stress 0.18963
## Run 77 stress 0.2034862
## Run 78 stress 0.2066994
## Run 79 stress 0.1954976
## Run 80 stress 0.2052672
## Run 81 stress 0.222485
## Run 82 stress 0.189617
## Run 83 stress 0.2036731
## Run 84 stress 0.2016689
## Run 85 stress 0.1884583
## ... Procrustes: rmse 0.02027165 max resid 0.1223164
## Run 86 stress 0.2166119
## Run 87 stress 0.2075466
## Run 88 stress 0.2079606
## Run 89 stress 0.1886734
## Run 90 stress 0.2018272
## Run 91 stress 0.1966538
## Run 92 stress 0.2178612
## Run 93 stress 0.1888607
## Run 94 stress 0.1977504
## Run 95 stress 0.2016763
## Run 96 stress 0.1881814
## ... Procrustes: rmse 0.0319324 max resid 0.1595439
## Run 97 stress 0.2066725
## Run 98 stress 0.2061703
## Run 99 stress 0.2073608
## Run 100 stress 0.1932321
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 97: stress ratio > sratmax
## 3: scale factor of the gradient < sfgrminanosim_no_trees <- anosim(wisconsin(sqrt(plant_com_no_trees)), plant_com_env$Type, distance = "bray", permutations = 9999)
# Create label for the plot
R_val <- round(anosim_no_trees$statistic, 3)
p_val <- ifelse(anosim_no_trees$signif < 0.001, "< 0.001", paste0("= ", signif(anosim_no_trees$signif, 2)))
anosim_label <- paste0("ANOSIM R = ", R_val, "\n", "p ", p_val)ordi_ellipses_no_trees <- ordiellipse(NMDS_plant_com_no_trees, plant_com_env$Type, kind = "se", conf = 0.95, draw = "none")
ellipse_df_no_trees <- do.call(rbind, lapply(names(ordi_ellipses_no_trees), function(group) {
ellipse_points <- vegan:::veganCovEllipse(ordi_ellipses_no_trees[[group]]$cov,
ordi_ellipses_no_trees[[group]]$center,
ordi_ellipses_no_trees[[group]]$scale)
data.frame(ellipse_points, Type = group)
}))
colnames(ellipse_df_no_trees)[1:2] <- c("NMDS1", "NMDS2")# Site scores
nmds_df_no_trees <- data.frame(
NMDS1 = NMDS_plant_com_no_trees$points[, 1],
NMDS2 = NMDS_plant_com_no_trees$points[, 2],
Type = plant_com_env$Type
)
# Species scores for labels (optional, may require handling NAs)
species_scores <- data.frame(scores(NMDS_plant_com_no_trees, display = "species"))
species_scores$Species <- rownames(species_scores)
colnames(species_scores)[1:2] <- c("NMDS1", "NMDS2")plant_nmds_no_trees_graph <- ggplot(nmds_df_no_trees, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
geom_point(shape = 21, size = 1, stroke = 1) +
geom_polygon(data = ellipse_df_no_trees, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
alpha = 0.4, linewidth = 0.5) +
geom_text(data = species_scores,
aes(x = NMDS1, y = NMDS2, label = Species),
inherit.aes = FALSE,
color = "white", size = 2, alpha = 0) +
scale_color_manual(name = "Forest Treatment",
values = colvec) +
scale_fill_manual(name = "Forest Treatment",
values = colvec) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10),) +
annotate("text", x = Inf, y = Inf, label = anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
labs(title = "Plant Community Separation (Without Tree Species)")
plant_nmds_no_trees_graph <- plant_nmds_no_trees_graph + coord_fixed()
plant_nmds_no_trees_graph
save_plot(plant_nmds_no_trees_graph, "plant_nmds_no_trees_graph")## png
## 2Plant NMDS resampling (No Trees)
# Base ANOSIM R value with no tree species
base_R_no_trees <- anosim(wisconsin(sqrt(plant_com_no_trees)), plant_com_env$Type, distance = "bray", permutations = 999)$statistic
# Store results
jackknife_results_no_trees <- data.frame(Species = character(),
R_value = numeric(),
Delta_R = numeric(),
stringsAsFactors = FALSE)
# Loop over all species (columns) in the no-trees matrix
species_names_no_trees <- colnames(plant_com_no_trees)
for (sp in species_names_no_trees) {
# Remove one species
reduced_matrix <- plant_com_no_trees[, setdiff(species_names_no_trees, sp)]
# Run ANOSIM
anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), plant_com_env$Type, distance = "bray", permutations = 999)
# Store R and change from base
R_val <- anosim_result$statistic
delta_R <- base_R_no_trees - R_val
jackknife_results_no_trees <- rbind(jackknife_results_no_trees,
data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}
# Sort by greatest impact (largest Delta_R)
jackknife_results_no_trees <- jackknife_results_no_trees[order(-jackknife_results_no_trees$Delta_R), ]
# Plot the results
plant_jk_no_trees_graph <- ggplot(jackknife_results_no_trees, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
geom_col(fill = "#0072B2") +
coord_flip() +
labs(title = "Plant Species Impact on NMDS Clustering (No Trees)",
x = "Species Removed",
y = "Decrease in ANOSIM R") +
theme_minimal()
jackknife_results_no_trees## Species R_value Delta_R
## 44 Silver.birch 0.3659864 4.875403e-02
## 24 Hylocomium.splendens 0.3755024 3.923802e-02
## 3 Blaeberry 0.3780339 3.670650e-02
## 17 Downy.birch 0.3845283 3.021212e-02
## 2 Bent.sp. 0.3984823 1.625814e-02
## 10 Cowberry 0.4054596 9.280855e-03
## 12 Cross.leaved.heather 0.4090427 5.697729e-03
## 20 Heath.bedstraw 0.4107332 4.007247e-03
## 5 Bog.myrtle 0.4108004 3.940059e-03
## 33 Pleurozia.purpurea 0.4112755 3.464949e-03
## 19 Hard.fern 0.4120703 2.670098e-03
## 4 Bog.asphodel 0.4124495 2.290969e-03
## 15 Devil.s.bit.scabious 0.4127284 2.012022e-03
## 39 Rhytidiadelphus.loreus 0.4133907 1.349746e-03
## 36 Polytrichum.sp. 0.4136990 1.041404e-03
## 47 Star.sedge 0.4137026 1.037805e-03
## 42 Scots.pine 0.4137062 1.034206e-03
## 1 Bell.heather 0.4137752 9.652185e-04
## 13 Crowberry 0.4137758 9.646187e-04
## 32 Moss.sp. 0.4139216 8.188461e-04
## 9 Cotton.grass.sp. 0.4139330 8.074482e-04
## 21 Heath.milkwort 0.4140314 7.090667e-04
## 16 Devil.s.bit.scabious.1 0.4140961 6.442789e-04
## 41 Round.leaved.sundew 0.4143157 4.247202e-04
## 38 Racomitrium.sp. 0.4144393 3.011434e-04
## 26 Ling.heather 0.4146858 5.458974e-05
## 52 Unknown.moss.sp..3 0.4147848 -4.439165e-05
## 7 Club.moss 0.4148028 -6.238827e-05
## 28 Marsh.lousewort 0.4149264 -1.859650e-04
## 54 Wood.anemone 0.4151759 -4.355181e-04
## 11 Cowslip 0.4152011 -4.607134e-04
## 56 Woodrush.sp. 0.4154225 -6.820718e-04
## 27 Liverwort.sp. 0.4157674 -1.027007e-03
## 25 Lemon.scented.fern 0.4158274 -1.086996e-03
## 8 Common.cow.wheat 0.4158874 -1.146984e-03
## 37 Ptilium.crista.castrensis 0.4162186 -1.478122e-03
## 51 Unknown.moss.sp..2 0.4162515 -1.511116e-03
## 18 Eyebright 0.4163085 -1.568105e-03
## 40 Ribwort.plantain 0.4165155 -1.775066e-03
## 23 Horsetail.sp. 0.4165377 -1.797262e-03
## 49 Tufted.hairgrass 0.4165377 -1.797262e-03
## 22 Heath.speedwell 0.4167656 -2.025219e-03
## 29 Meadow.buttercup 0.4168730 -2.132599e-03
## 35 Poa.sp. 0.4169684 -2.227981e-03
## 55 Wood.sorrel 0.4170224 -2.281971e-03
## 30 Meadow.grass.sp. 0.4170872 -2.346759e-03
## 6 Bracken 0.4172384 -2.497930e-03
## 14 Deergrass 0.4175659 -2.825469e-03
## 50 Unknown.moss.sp. 0.4178514 -3.111015e-03
## 45 Soft.rush 0.4180110 -3.270585e-03
## 43 Sedge.sp. 0.4181508 -3.410359e-03
## 53 Violet.sp. 0.4189186 -4.178214e-03
## 34 Pleurozium.schreberi 0.4198161 -5.075646e-03
## 31 Molinia.caerulea 0.4229829 -8.242450e-03
## 48 Tormentil 0.4470726 -3.233212e-02
## 46 Sphagnum.sp. 0.4488872 -3.414678e-02# Display the graph
plant_jk_no_trees_graph
save_plot(plant_jk_no_trees_graph, "plant_jk_no_trees")## png
## 2palette.colors()## [1] "#000000" "#E69F00" "#56B4E9" "#009E73" "#F0E442" "#0072B2" "#D55E00"
## [8] "#CC79A7" "#999999"Soil Invertebrate Diversity:
Extracting model predicted data
predicted_soil_inv_diversity <- predict(soil_invert_diversity_model, type = "response", level = 1)
predicted_soil_div_with_ci <- intervals(soil_invert_diversity_model, level = 0.95)
# Extract just the fixed effects into a data frame
soil_div_ci_df <- as.data.frame(predicted_soil_div_with_ci$fixed)
# Optionally, name the columns
colnames(soil_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")
# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
soil_div_ci_df$Effect <- rownames(predicted_soil_div_with_ci$fixed)
# Reorder for readability
soil_div_ci_df <- soil_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]
# View it
print(soil_div_ci_df)## Effect Estimate CI_lower
## (Intercept) (Intercept) 0.36597967 0.26761870
## TypeRegen TypeRegen 0.09497504 -0.04674363
## TypeUnforested (heath) TypeUnforested (heath) -0.06544998 -0.21779287
## TypeUnforested (remnant) TypeUnforested (remnant) -0.05195028 -0.19366895
## CI_upper
## (Intercept) 0.46434064
## TypeRegen 0.23669371
## TypeUnforested (heath) 0.08689291
## TypeUnforested (remnant) 0.08976839# Extract the intercept value
soil_intercept_estimate <- soil_div_ci_df$Estimate[soil_div_ci_df$Effect == "(Intercept)"]
soil_intercept_lower <- soil_div_ci_df$CI_lower[soil_div_ci_df$Effect == "(Intercept)"]
soil_intercept_upper <- soil_div_ci_df$CI_upper[soil_div_ci_df$Effect == "(Intercept)"]
# Now create a new column with actual predicted values
soil_div_ci_df$Actual_Estimate <- soil_div_ci_df$Estimate
soil_div_ci_df$Actual_CI_lower <- soil_div_ci_df$CI_lower
soil_div_ci_df$Actual_CI_upper <- soil_div_ci_df$CI_upper
# Add the intercept to all *non-intercept* rows
soil_non_intercept_rows <- soil_div_ci_df$Effect != "(Intercept)"
soil_div_ci_df$Actual_Estimate[soil_non_intercept_rows] <-
soil_div_ci_df$Estimate[soil_non_intercept_rows] + soil_intercept_estimate
soil_div_ci_df$Actual_CI_lower[soil_non_intercept_rows] <-
soil_div_ci_df$CI_lower[soil_non_intercept_rows] + soil_intercept_lower
soil_div_ci_df$Actual_CI_upper[soil_non_intercept_rows] <-
soil_div_ci_df$CI_upper[soil_non_intercept_rows] + soil_intercept_upper
# Using base R
colnames(soil_div_ci_df)[colnames(soil_div_ci_df) == "Effect"] <- "Type"
soil_div_ci_df$Type[soil_div_ci_df$Type == "(Intercept)"] <- "Mature"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeRegen"] <- "Regen"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
soil_div_ci_df## Type Estimate CI_lower
## (Intercept) Mature 0.36597967 0.26761870
## TypeRegen Regen 0.09497504 -0.04674363
## TypeUnforested (heath) Unforested (heath) -0.06544998 -0.21779287
## TypeUnforested (remnant) Unforested (remnant) -0.05195028 -0.19366895
## CI_upper Actual_Estimate Actual_CI_lower
## (Intercept) 0.46434064 0.3659797 0.26761870
## TypeRegen 0.23669371 0.4609547 0.22087506
## TypeUnforested (heath) 0.08689291 0.3005297 0.04982583
## TypeUnforested (remnant) 0.08976839 0.3140294 0.07394974
## Actual_CI_upper
## (Intercept) 0.4643406
## TypeRegen 0.7010343
## TypeUnforested (heath) 0.5512335
## TypeUnforested (remnant) 0.5541090Plotting soil invertebrate diversity
soil_diversity_graph <- ggplot(soil_invert_div, aes(x = Type, y = Diversity_Index, fill = Type)) +
# Violin plot to show distribution of raw diversity index values
geom_violin(alpha = 0.3, color = NA, width = 0.9) +
# Individual data points (jittered)
geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
# Model-estimated means (points)
geom_point(data = soil_div_ci_df,
aes(x = Type, y = Actual_Estimate),
shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +
# Confidence intervals (error bars)
geom_errorbar(data = soil_div_ci_df,
aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper),
width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +
# Theme and styling
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10),
legend.position = "none"
) +
labs(
x = "Forest Treatment",
y = "Diversity Index (1 - Simpson's D)",
fill = "Type",
title = "Soil Invertebrate Diversity"
) +
scale_fill_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))Soil Invertebrate Richness:
Extracting soil invertebrate richness data:
# Step 1: Extract the intercept value
soil_intercept_value <- coef(summary(soil_invert_richness_model))["(Intercept)", "Estimate"]
soil_wald_ci <- confint(soil_invert_richness_model, level = 0.95)## Waiting for profiling to be done...# Drop .sig01 row
soil_wald_ci <- soil_wald_ci[!rownames(soil_wald_ci) %in% ".sig01", ]
# Then build the dataframe
soil_results_df <- data.frame(
Effect = rownames(coef(summary(soil_invert_richness_model))),
Estimate = coef(summary(soil_invert_richness_model))[, "Estimate"],
Std_Error = coef(summary(soil_invert_richness_model))[, "Std. Error"],
p_value = coef(summary(soil_invert_richness_model))[, "Pr(>|z|)"],
CI_Lower = soil_wald_ci[, 1],
CI_Upper = soil_wald_ci[, 2]
)
# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
soil_results_df$Real_Estimate <- soil_results_df$Estimate
soil_results_df$Real_CI_Lower <- soil_results_df$CI_Lower
soil_results_df$Real_CI_Upper <- soil_results_df$CI_Upper
soil_non_intercepts <- soil_results_df$Effect != "(Intercept)"
soil_results_df$Real_Estimate[soil_non_intercepts] <- soil_results_df$Estimate[soil_non_intercepts] + soil_intercept_value
soil_results_df$Real_CI_Lower[soil_non_intercepts] <- soil_results_df$CI_Lower[soil_non_intercepts] + soil_intercept_value
soil_results_df$Real_CI_Upper[soil_non_intercepts] <- soil_results_df$CI_Upper[soil_non_intercepts] + soil_intercept_value
# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(soil_results_df)[colnames(soil_results_df) == "Effect"] <- "Type"
soil_results_df$Type[soil_results_df$Type == "(Intercept)"] <- "Mature"
soil_results_df$Type[soil_results_df$Type == "TypeRegen"] <- "Regen"
soil_results_df$Type[soil_results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
soil_results_df$Type[soil_results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(soil_results_df) <- soil_results_df$Type
# Step 8: Exponentiate Estimates and Confidence Intervals
soil_results_df$Real_Estimate <- exp(soil_results_df$Real_Estimate)
soil_results_df$Real_CI_Lower <- exp(soil_results_df$Real_CI_Lower)
soil_results_df$Real_CI_Upper <- exp(soil_results_df$Real_CI_Upper)
print(soil_results_df)## Type Estimate Std_Error p_value
## Mature Mature 1.3312346 0.1373606 3.275698e-22
## Regen Regen -0.2852660 0.2142310 1.829978e-01
## Unforested (heath) Unforested (heath) -0.3379828 0.2364418 1.528736e-01
## Unforested (remnant) Unforested (remnant) -0.3996764 0.2217452 7.148053e-02
## CI_Lower CI_Upper Real_Estimate Real_CI_Lower
## Mature 1.0493695 1.58889621 3.785714 2.855850
## Regen -0.7119726 0.13093157 2.846154 1.857557
## Unforested (heath) -0.8152634 0.11588244 2.700000 1.675265
## Unforested (remnant) -0.8434047 0.02932028 2.538462 1.628778
## Real_CI_Upper
## Mature 4.898339
## Regen 4.315297
## Unforested (heath) 4.250842
## Unforested (remnant) 3.898356Plotting soil invertebrate richness
soil_invert_richness_graph <- ggplot(soil_invert_rich, aes(x = Type, y = Richness, fill = Type)) +
# Violin plot for distribution of raw richness values
geom_violin(alpha = 0.3, color = NA, width = 0.9) +
# Jittered individual data points
geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
# Model-based means (points)
geom_point(data = soil_results_df,
aes(x = Type, y = Real_Estimate),
shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +
# Model-based confidence intervals (error bars)
geom_errorbar(data = soil_results_df,
aes(x = Type, ymin = Real_CI_Lower, ymax = Real_CI_Upper),
width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +
# Clean minimal theme
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10)
) +
labs(
x = "Forest Treatment",
y = "Order level Richness",
fill = "Type",
title = "Soil invertebrate Richness"
) +
scale_fill_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))
soil_div_rich_graph <- soil_diversity_graph + soil_invert_richness_graph
soil_div_rich_graph
save_plot(soil_div_rich_graph, "soil_rich_div_graph")## png
## 2Soil Invertebrate Community Composition:
#soil_invert_dataset <- plant_com_dataset[-c(561, 701), ]
soil_invert_dataset <- soil_invert_dataset[c("Order", "Type", "Abundance", "Site.y", "Plot")]
soil_invert_dataset <- aggregate(cbind(Abundance) ~ Order + Plot, data = soil_invert_dataset, FUN = sum)
soil_invert_dataset <- merge(soil_invert_dataset, overview_dataset, by = "Plot")
soil_invert_dataset <- soil_invert_dataset[c("Order", "Type", "Abundance", "Site", "Plot")]
soil_invert_wide <- as.data.frame(spread(soil_invert_dataset, Order, Abundance))After converting the data into wide format, I need to replace all the NA values with 0
soil_invert_wide[is.na(soil_invert_wide)] <- 0
soil_invert_env <- soil_invert_wide[c("Plot", "Type", "Site")]
soil_invert_wide$Plot <- NULL
soil_invert_wide$Site <- NULL
soil_invert_wide$Type <- NULL
soil_invert_wide$Family <- NULL
soil_invert_wide$Genus <- NULL
soil_invert_wide <- as.data.frame(lapply(soil_invert_wide, as.integer))
str(soil_invert_wide)## 'data.frame': 50 obs. of 11 variables:
## $ Acari : int 11 14 217 75 304 10 55 70 109 27 ...
## $ Araneae : int 0 0 0 0 0 0 0 0 0 1 ...
## $ Coleoptera : int 0 0 2 0 0 0 0 0 0 0 ...
## $ Diptera : int 4 1 2 1 4 5 4 2 1 8 ...
## $ Entomobryomorpha: int 0 0 0 0 0 1 0 0 0 0 ...
## $ Geophilomorpha : int 0 0 0 0 0 0 0 1 0 0 ...
## $ Hemiptera : int 0 0 0 0 10 0 0 0 0 0 ...
## $ Hymenoptera : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Isopoda : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Podomorpha : int 5 1 70 9 46 0 7 61 9 3 ...
## $ Thysanoptera : int 0 0 0 0 0 0 0 0 0 0 ...NMDS_soil_invert <- metaMDS(soil_invert_wide)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1761955
## Run 1 stress 0.1712266
## ... New best solution
## ... Procrustes: rmse 0.04440726 max resid 0.2888872
## Run 2 stress 0.2220741
## Run 3 stress 0.1926232
## Run 4 stress 0.2231124
## Run 5 stress 0.1725766
## Run 6 stress 0.1836321
## Run 7 stress 0.1843298
## Run 8 stress 0.176317
## Run 9 stress 0.1730456
## Run 10 stress 0.197426
## Run 11 stress 0.1884122
## Run 12 stress 0.2083265
## Run 13 stress 0.1992799
## Run 14 stress 0.1906248
## Run 15 stress 0.2098807
## Run 16 stress 0.1727093
## Run 17 stress 0.1744323
## Run 18 stress 0.1773268
## Run 19 stress 0.1758481
## Run 20 stress 0.1776654
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmaxNMDS_soil_invert##
## Call:
## metaMDS(comm = soil_invert_wide)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(soil_invert_wide))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1712266
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 1 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(soil_invert_wide))'stressplot(NMDS_soil_invert)
soil_invert_env## Plot Type Site
## 1 63 Mature Glen Cluanie
## 2 66 Mature Glen Cluanie
## 3 68 Mature Glen Cluanie
## 4 69 Mature Glen Cluanie
## 5 74 Mature Glen Cluanie
## 6 81 Mature Glen Cluanie
## 7 82 Mature Glen Cluanie
## 8 3 Mature Glen Mallie
## 9 4 Mature Glen Mallie
## 10 10 Mature Glen Mallie
## 11 12 Mature Glen Mallie
## 12 14 Mature Glen Mallie
## 13 24 Mature Glen Mallie
## 14 25 Mature Glen Mallie
## 15 62 Regen Glen Cluanie
## 16 64 Regen Glen Cluanie
## 17 65 Regen Glen Cluanie
## 18 67 Regen Glen Cluanie
## 19 71 Regen Glen Cluanie
## 20 80 Regen Glen Cluanie
## 21 6 Regen Glen Mallie
## 22 7 Regen Glen Mallie
## 23 8 Regen Glen Mallie
## 24 11 Regen Glen Mallie
## 25 13 Regen Glen Mallie
## 26 23 Regen Glen Mallie
## 27 26 Regen Glen Mallie
## 28 57 Unforested (heath) Glen Cluanie
## 29 58 Unforested (heath) Glen Cluanie
## 30 59 Unforested (heath) Glen Cluanie
## 31 60 Unforested (heath) Glen Cluanie
## 32 77 Unforested (heath) Glen Cluanie
## 33 16 Unforested (heath) Glen Mallie
## 34 17 Unforested (heath) Glen Mallie
## 35 18 Unforested (heath) Glen Mallie
## 36 19 Unforested (heath) Glen Mallie
## 37 20 Unforested (heath) Glen Mallie
## 38 61 Unforested (remnant) Glen Cluanie
## 39 70 Unforested (remnant) Glen Cluanie
## 40 72 Unforested (remnant) Glen Cluanie
## 41 75 Unforested (remnant) Glen Cluanie
## 42 76 Unforested (remnant) Glen Cluanie
## 43 78 Unforested (remnant) Glen Cluanie
## 44 1 Unforested (remnant) Glen Mallie
## 45 2 Unforested (remnant) Glen Mallie
## 46 5 Unforested (remnant) Glen Mallie
## 47 9 Unforested (remnant) Glen Mallie
## 48 15 Unforested (remnant) Glen Mallie
## 49 22 Unforested (remnant) Glen Mallie
## 50 27 Unforested (remnant) Glen Malliesoil_invert_env$Type <- factor(soil_invert_env$Type)
scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_soil_invert, type = "n")
with(soil_invert_env, points(NMDS_soil_invert, display = "sites", col = colvec[Type],
scaling = scl, pch = 21, bg = colvec[Type]))## Warning in plot.xy(xy.coords(x, y), type = type, ...): "scaling" is not a
## graphical parameterwith(soil_invert_env, legend("topleft", legend = levels(Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec))
NMDS_soil_invert.fit <- envfit(NMDS_soil_invert ~ Type , data = soil_invert_env, perm = 999)
NMDS_soil_invert.fit##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## TypeMature -0.0571 0.2111
## TypeRegen -0.0450 0.1406
## TypeUnforested (heath) -0.0242 -0.1613
## TypeUnforested (remnant) 0.1251 -0.2439
##
## Goodness of fit:
## r2 Pr(>r)
## Type 0.0823 0.224
## Permutation: free
## Number of permutations: 999plot(NMDS_soil_invert.fit, dis = "site")## Warning in text.default(x$factors$centroids[, choices, drop = FALSE], labs$f, :
## "dis" is not a graphical parameterplot(NMDS_soil_invert.fit)
### Soil invertebrate NMDS
# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)
soil_ordi_ellipses <- ordiellipse(NMDS_soil_invert, soil_invert_env$Type, kind = "se", conf = 0.95, draw = "none")
# 2. Extract ellipse coordinates into a data frame
soil_ellipse_df <- do.call(rbind, lapply(names(soil_ordi_ellipses), function(group) {
soil_ellipse_points <- vegan:::veganCovEllipse(soil_ordi_ellipses[[group]]$cov,
soil_ordi_ellipses[[group]]$center,
soil_ordi_ellipses[[group]]$scale)
data.frame(soil_ellipse_points, Type = group)
}))
colnames(soil_ellipse_df)[1:2] <- c("NMDS1", "NMDS2")
library(ggplot2)
# Your main NMDS points
soil_nmds_df <- data.frame(
NMDS1 = NMDS_soil_invert$points[, 1],
NMDS2 = NMDS_soil_invert$points[, 2],
Type = soil_invert_env$Type
)
soil_anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(soil_invert_wide)), soil_invert_env$Type, distance = "bray", permutations = 9999)
soil_R_val <- round(soil_anosim_plant_com_r_tr$statistic, 3)
soil_p_val <- ifelse(soil_anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(soil_anosim_plant_com_r_tr$signif, 2)))
soil_anosim_label <- paste0("ANOSIM R = ", soil_R_val, "\n", "p ", soil_p_val)
# Site scores
soil_site_scores <- as.data.frame(NMDS_soil_invert$points)
colnames(soil_site_scores)[1:2] <- c("NMDS1", "NMDS2")
soil_site_scores$Type <- soil_invert_env$Type # Add grouping
# Species scores (usually in $species, or use scores())
soil_species_scores <- as.data.frame(scores(NMDS_soil_invert, display = "species"))
soil_species_scores$Species <- rownames(soil_species_scores)
soil_nmds_hull_graph <- ggplot(soil_nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
geom_point(shape = 21, size = 1, stroke = 1) +
geom_polygon(data = soil_ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
alpha = 0.4, linewidth = 0.5) +
# Species labels
geom_text(data = soil_species_scores, aes(x = NMDS1, y = NMDS2, label = Species),
color = "darkcyan", alpha = 0,
inherit.aes = FALSE,
size = 3) +
scale_color_manual(name = "Forest Treatment",
values = colvec) +
scale_fill_manual(name = "Forest Treatment",
values = colvec) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10)) +
annotate("text", x = Inf, y = Inf, label = soil_anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
labs(title = "Soil Invertebrate Community Separation")
soil_nmds_hull_graph + coord_fixed()
save_plot(soil_nmds_hull_graph, "soil_nmds_hull_graph")## png
## 2Soil invertebrate NMDS resampling
library(vegan)
soil_invert_wide_nonzero <- soil_invert_wide[rowSums(soil_invert_wide) > 0, ]
soil_invert_env_nonzero <- soil_invert_env[rowSums(soil_invert_wide) > 0, ]
nonzero_rows <- rowSums(soil_invert_wide) > 0
soil_invert_wide_nonzero <- soil_invert_wide[nonzero_rows, ]
soil_invert_env_nonzero <- soil_invert_env[nonzero_rows, ]
anyNA(soil_invert_wide_nonzero) # Should return FALSE## [1] FALSElibrary(vegan)
base_R <- anosim(
wisconsin(sqrt(soil_invert_wide_nonzero)),
soil_invert_env_nonzero$Type,
distance = "bray",
permutations = 999
)$statistic
# Use non-zero rows only
base_R <- anosim(wisconsin(sqrt(soil_invert_wide_nonzero)), soil_invert_env_nonzero$Type, distance = "bray", permutations = 999)$statistic
jackknife_results <- data.frame(Species = character(),
R_value = numeric(),
Delta_R = numeric(),
stringsAsFactors = FALSE)
species_names <- colnames(soil_invert_wide_nonzero)
for (sp in species_names) {
# Remove the species (column)
reduced_matrix <- soil_invert_wide_nonzero[, setdiff(species_names, sp)]
# Remove rows that are now completely zero
nonzero_rows <- rowSums(reduced_matrix) > 0
reduced_matrix <- reduced_matrix[nonzero_rows, ]
reduced_env <- soil_invert_env_nonzero[nonzero_rows, ]
# Skip iteration if data becomes too small or empty
if (nrow(reduced_matrix) < 2 || anyNA(reduced_matrix)) next
# Run ANOSIM
anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), reduced_env$Type, distance = "bray", permutations = 999)
# Store results
R_val <- anosim_result$statistic
delta_R <- base_R - R_val
jackknife_results <- rbind(jackknife_results, data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}
# Sort
jackknife_results <- jackknife_results[order(-jackknife_results$Delta_R), ]
# Plot
soil_inv_jk_graph <- ggplot(jackknife_results, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
geom_col(fill = "#0072B2") +
coord_flip() +
labs(title = "Soil Invertebrate Order Impact on NMDS Clustering",
x = "Order Removed",
y = "Decrease in ANOSIM R") +
theme_minimal()
soil_inv_jk_graph
save_plot(soil_inv_jk_graph, "soil_inv_jk_graph")## png
## 2Ground Invertebrate Diversity:
Extracting model predicted ground invertebrate diversity:
predicted_pitfall_diversity <- predict(pitfall_diversity_model, type = "response", level = 1)
predicted_pitfall_div_with_ci <- intervals(pitfall_diversity_model, level = 0.95)
# Extract just the fixed effects into a data frame
pit_div_ci_df <- as.data.frame(predicted_pitfall_div_with_ci$fixed)
# Optionally, name the columns
colnames(pit_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")
# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
pit_div_ci_df$Effect <- rownames(predicted_pitfall_div_with_ci$fixed)
# Reorder for readability
pit_div_ci_df <- pit_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]
# View it
print(pit_div_ci_df)## Effect Estimate CI_lower
## (Intercept) (Intercept) 0.661053722 0.55147425
## TypeRegen TypeRegen -0.002566658 -0.12499009
## TypeUnforested (heath) TypeUnforested (heath) 0.077685011 -0.05544383
## TypeUnforested (remnant) TypeUnforested (remnant) -0.004524457 -0.12391595
## CI_upper
## (Intercept) 0.7706332
## TypeRegen 0.1198568
## TypeUnforested (heath) 0.2108138
## TypeUnforested (remnant) 0.1148670# Extract the intercept value
pit_intercept_estimate <- pit_div_ci_df$Estimate[pit_div_ci_df$Effect == "(Intercept)"]
pit_intercept_lower <- pit_div_ci_df$CI_lower[pit_div_ci_df$Effect == "(Intercept)"]
pit_intercept_upper <- pit_div_ci_df$CI_upper[pit_div_ci_df$Effect == "(Intercept)"]
# Now create a new column with actual predicted values
pit_div_ci_df$Actual_Estimate <- pit_div_ci_df$Estimate
pit_div_ci_df$Actual_CI_lower <- pit_div_ci_df$CI_lower
pit_div_ci_df$Actual_CI_upper <- pit_div_ci_df$CI_upper
# Add the intercept to all *non-intercept* rows
pit_non_intercept_rows <- pit_div_ci_df$Effect != "(Intercept)"
pit_div_ci_df$Actual_Estimate[pit_non_intercept_rows] <-
pit_div_ci_df$Estimate[pit_non_intercept_rows] + pit_intercept_estimate
pit_div_ci_df$Actual_CI_lower[pit_non_intercept_rows] <-
pit_div_ci_df$CI_lower[pit_non_intercept_rows] + pit_intercept_lower
pit_div_ci_df$Actual_CI_upper[pit_non_intercept_rows] <-
pit_div_ci_df$CI_upper[pit_non_intercept_rows] + pit_intercept_upper
# Using base R
colnames(pit_div_ci_df)[colnames(pit_div_ci_df) == "Effect"] <- "Type"
pit_div_ci_df$Type[pit_div_ci_df$Type == "(Intercept)"] <- "Mature"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeRegen"] <- "Regen"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
pit_div_ci_df## Type Estimate CI_lower
## (Intercept) Mature 0.661053722 0.55147425
## TypeRegen Regen -0.002566658 -0.12499009
## TypeUnforested (heath) Unforested (heath) 0.077685011 -0.05544383
## TypeUnforested (remnant) Unforested (remnant) -0.004524457 -0.12391595
## CI_upper Actual_Estimate Actual_CI_lower
## (Intercept) 0.7706332 0.6610537 0.5514742
## TypeRegen 0.1198568 0.6584871 0.4264842
## TypeUnforested (heath) 0.2108138 0.7387387 0.4960304
## TypeUnforested (remnant) 0.1148670 0.6565293 0.4275583
## Actual_CI_upper
## (Intercept) 0.7706332
## TypeRegen 0.8904900
## TypeUnforested (heath) 0.9814470
## TypeUnforested (remnant) 0.8855002Plotting ground invertebrate diversity
pitfall_diversity_graph <- ggplot(pitfall_div, aes(x = Type, y = Diversity_Index, fill = Type)) +
# Violin plot for raw data distribution
geom_violin(alpha = 0.3, color = NA, width = 0.9) +
# Jittered individual data points
geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
# Model-estimated means
geom_point(data = pit_div_ci_df,
aes(x = Type, y = Actual_Estimate),
shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +
# Model-estimated confidence intervals
geom_errorbar(data = pit_div_ci_df,
aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper),
width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +
# Clean, minimal theme
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey", color = "black", size = 1),
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"),
panel.background = element_rect(fill = "white", color = "black"),
plot.margin = margin(10, 10, 10, 10),
legend.position = "none"
) +
labs(
x = "Forest Treatment",
y = "Diversity Index (1 - Simpson's D)",
fill = "Type",
title = "Ground Invertebrate Diversity"
) +
scale_fill_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))Pitfall Invertebrate Richness:
Extracting model predicted pitfall invertebrate richness:
# Step 1: Extract the intercept value
pit_intercept_value <- coef(summary(pitfall_richness_model))["(Intercept)", "Estimate"]
pit_wald_ci <- confint(pitfall_richness_model, level = 0.95)## Computing profile confidence intervals ...# Drop .sig01 row
pit_wald_ci <- pit_wald_ci[!rownames(pit_wald_ci) %in% ".sig01", ]
# Then build the dataframe
pit_results_df <- data.frame(
Effect = rownames(coef(summary(pitfall_richness_model))),
Estimate = coef(summary(pitfall_richness_model))[, "Estimate"],
Std_Error = coef(summary(pitfall_richness_model))[, "Std. Error"],
p_value = coef(summary(pitfall_richness_model))[, "Pr(>|z|)"],
CI_Lower = pit_wald_ci[, 1],
CI_Upper = pit_wald_ci[, 2]
)
# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
pit_results_df$Real_Estimate <- pit_results_df$Estimate
pit_results_df$Real_CI_Lower <- pit_results_df$CI_Lower
pit_results_df$Real_CI_Upper <- pit_results_df$CI_Upper
pit_non_intercepts <- pit_results_df$Effect != "(Intercept)"
pit_results_df$Real_Estimate[pit_non_intercepts] <- pit_results_df$Estimate[pit_non_intercepts] + pit_intercept_value
pit_results_df$Real_CI_Lower[pit_non_intercepts] <- pit_results_df$CI_Lower[pit_non_intercepts] + pit_intercept_value
pit_results_df$Real_CI_Upper[pit_non_intercepts] <- pit_results_df$CI_Upper[pit_non_intercepts] + pit_intercept_value
# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(pit_results_df)[colnames(pit_results_df) == "Effect"] <- "Type"
pit_results_df$Type[pit_results_df$Type == "(Intercept)"] <- "Mature"
pit_results_df$Type[pit_results_df$Type == "TypeRegen"] <- "Regen"
pit_results_df$Type[pit_results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
pit_results_df$Type[pit_results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"
# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(pit_results_df) <- pit_results_df$Type
# Step 8: Exponentiate Estimates and Confidence Intervals
pit_results_df$Real_Estimate <- exp(pit_results_df$Real_Estimate)
pit_results_df$Real_CI_Lower <- exp(pit_results_df$Real_CI_Lower)
pit_results_df$Real_CI_Upper <- exp(pit_results_df$Real_CI_Upper)
print(pit_results_df)## Type Estimate Std_Error p_value
## Mature Mature 1.7717679 0.1278215 1.087372e-43
## Regen Regen -0.2433485 0.1347723 7.097659e-02
## Unforested (heath) Unforested (heath) -0.1389162 0.1539754 3.669519e-01
## Unforested (remnant) Unforested (remnant) -0.1734879 0.1337063 1.944492e-01
## CI_Lower CI_Upper Real_Estimate Real_CI_Lower
## Mature 1.4471914 2.08273351 5.881242 4.251158
## Regen -0.5095904 0.02007054 4.610883 3.533106
## Unforested (heath) -0.4462172 0.15916599 5.118450 3.764258
## Unforested (remnant) -0.4374303 0.08802427 4.944521 3.797480
## Real_CI_Upper
## Mature 8.026379
## Regen 6.000474
## Unforested (heath) 6.895947
## Unforested (remnant) 6.422402Plotting ground invertebrate richness
pitfall_richness_graph <- ggplot(pit_results_df, aes(x = Type, y = Real_Estimate, fill = Type)) +
geom_jitter(data = pitfall_rich, aes(x = Type, y = Richness, color = Type),
size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) + # Add individual data points
geom_point(size = 3, shape = 21, color = "black", fill = "black", position = position_dodge(width = 0.8)) + # Scatter plot with mean
geom_errorbar(
aes(ymin = Real_CI_Lower, ymax = Real_CI_Upper), # Error bars based on confidence intervals
width = 0.2, # Error bar width
color = "black", # Error bar color
size = 0.5 # Error bar size
) +
geom_violin(data = pitfall_rich, aes(x = Type, y = Richness, fill = Type),
alpha = 0.3, color = NA, width = 0.9) +
theme_bw() + # Minimal theme
theme(
axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
strip.background = element_rect(fill = "grey", color = "black", size = 1), # Grey background for facet labels
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"), # Add spacing between panels if needed
panel.background = element_rect(fill = "white", color = "black"), # Apply white background with black border to panels
plot.margin = margin(10, 10, 10, 10) # Text color inside the grey box
) +
labs(
x = "Forest Treatment", # x-axis label
y = "Order level Richness", # y-axis label
fill = "Type",
title = "Ground Invertebrate Richness" # Overall title
) +
scale_fill_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
)) +
scale_color_manual(
name = "Forest Treatment",
values = c(
"Mature" = "#D55E00",
"Regen" = "#009E73",
"Unforested (heath)" = "#56B4E9",
"Unforested (remnant)" = "#CC79A7"
))
pitfall_rich_div_graph <- pitfall_diversity_graph + pitfall_richness_graph
pitfall_rich_div_graph
save_plot(pitfall_rich_div_graph, "pitfall_rich_div_graph")## png
## 2Pitfall Community Composition:
pitfall_dataset <- pitfall_dataset[c("Order", "Type", "Abundance", "Site.y", "Plot")]
pitfall_dataset <- aggregate(cbind(Abundance) ~ Order + Plot, data = pitfall_dataset, FUN = sum)
pitfall_dataset <- merge(pitfall_dataset, overview_dataset, by = "Plot")
pitfall_dataset <- pitfall_dataset[c("Order", "Type", "Abundance", "Site", "Plot")]
pitfall_wide <- as.data.frame(spread(pitfall_dataset, Order, Abundance))After converting the data into wide format, I need to replace all the NA values with 0
pitfall_wide[is.na(pitfall_wide)] <- 0
pitfall_env <- pitfall_wide[c("Plot", "Type", "Site")]
pitfall_wide$Plot <- NULL
pitfall_wide$Site <- NULL
pitfall_wide$Type <- NULL
pitfall_wide$Family <- NULL
pitfall_wide$Genus <- NULL
pitfall_wide <- as.data.frame(lapply(pitfall_wide, as.integer))
str(pitfall_wide)## 'data.frame': 74 obs. of 14 variables:
## $ Acari : int 0 0 1 1 0 0 1 0 1 3 ...
## $ Araneae : int 0 2 3 1 0 1 0 2 3 1 ...
## $ Arionidae : int 0 0 0 1 0 0 0 0 0 0 ...
## $ Coleoptera : int 0 0 2 6 13 1 3 1 7 6 ...
## $ Diptera : int 2 0 2 5 1 10 2 10 2 5 ...
## $ Entomobryomorpha: int 2 0 3 1 2 7 1 6 5 14 ...
## $ Hemiptera : int 0 0 0 3 0 0 0 0 1 1 ...
## $ Hymenoptera : int 0 2 1 0 1 2 0 3 0 6 ...
## $ Isopoda : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Julida : int 0 1 0 1 0 0 0 0 0 0 ...
## $ Lepidoptera : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Opisthopora : int 0 0 1 1 0 0 1 0 0 0 ...
## $ Podomorpha : int 0 0 0 0 0 5 1 0 0 1 ...
## $ Symphpleona : int 0 0 0 0 0 0 0 0 0 0 ...NMDS_pitfall <- metaMDS(pitfall_wide)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2371
## Run 1 stress 0.2403582
## Run 2 stress 0.2379947
## Run 3 stress 0.2325491
## ... New best solution
## ... Procrustes: rmse 0.08307471 max resid 0.3866432
## Run 4 stress 0.2327488
## ... Procrustes: rmse 0.01792176 max resid 0.1000026
## Run 5 stress 0.242124
## Run 6 stress 0.2391028
## Run 7 stress 0.2325653
## ... Procrustes: rmse 0.01496153 max resid 0.1016296
## Run 8 stress 0.2329172
## ... Procrustes: rmse 0.0300439 max resid 0.1584433
## Run 9 stress 0.2331662
## Run 10 stress 0.2468438
## Run 11 stress 0.233911
## Run 12 stress 0.238065
## Run 13 stress 0.2327124
## ... Procrustes: rmse 0.04351231 max resid 0.2047893
## Run 14 stress 0.2417877
## Run 15 stress 0.2331679
## Run 16 stress 0.2326182
## ... Procrustes: rmse 0.01523301 max resid 0.1015965
## Run 17 stress 0.2349177
## Run 18 stress 0.2331387
## Run 19 stress 0.2321276
## ... New best solution
## ... Procrustes: rmse 0.02328379 max resid 0.103097
## Run 20 stress 0.2624901
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 3: no. of iterations >= maxit
## 17: stress ratio > sratmaxNMDS_pitfall##
## Call:
## metaMDS(comm = pitfall_wide)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(pitfall_wide))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2321276
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 19 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(pitfall_wide))'stressplot(NMDS_pitfall)
pitfall_env## Plot Type Site
## 1 63 Mature Glen Cluanie
## 2 66 Mature Glen Cluanie
## 3 68 Mature Glen Cluanie
## 4 69 Mature Glen Cluanie
## 5 74 Mature Glen Cluanie
## 6 81 Mature Glen Cluanie
## 7 82 Mature Glen Cluanie
## 8 36 Mature Glen Loyne
## 9 39 Mature Glen Loyne
## 10 40 Mature Glen Loyne
## 11 41 Mature Glen Loyne
## 12 43 Mature Glen Loyne
## 13 48 Mature Glen Loyne
## 14 50 Mature Glen Loyne
## 15 3 Mature Glen Mallie
## 16 4 Mature Glen Mallie
## 17 10 Mature Glen Mallie
## 18 12 Mature Glen Mallie
## 19 14 Mature Glen Mallie
## 20 24 Mature Glen Mallie
## 21 25 Mature Glen Mallie
## 22 62 Regen Glen Cluanie
## 23 64 Regen Glen Cluanie
## 24 65 Regen Glen Cluanie
## 25 67 Regen Glen Cluanie
## 26 71 Regen Glen Cluanie
## 27 73 Regen Glen Cluanie
## 28 80 Regen Glen Cluanie
## 29 35 Regen Glen Loyne
## 30 37 Regen Glen Loyne
## 31 38 Regen Glen Loyne
## 32 42 Regen Glen Loyne
## 33 44 Regen Glen Loyne
## 34 45 Regen Glen Loyne
## 35 47 Regen Glen Loyne
## 36 6 Regen Glen Mallie
## 37 7 Regen Glen Mallie
## 38 8 Regen Glen Mallie
## 39 11 Regen Glen Mallie
## 40 13 Regen Glen Mallie
## 41 23 Regen Glen Mallie
## 42 26 Regen Glen Mallie
## 43 57 Unforested (heath) Glen Cluanie
## 44 58 Unforested (heath) Glen Cluanie
## 45 59 Unforested (heath) Glen Cluanie
## 46 60 Unforested (heath) Glen Cluanie
## 47 77 Unforested (heath) Glen Cluanie
## 48 29 Unforested (heath) Glen Loyne
## 49 30 Unforested (heath) Glen Loyne
## 50 31 Unforested (heath) Glen Loyne
## 51 53 Unforested (heath) Glen Loyne
## 52 54 Unforested (heath) Glen Loyne
## 53 16 Unforested (heath) Glen Mallie
## 54 20 Unforested (heath) Glen Mallie
## 55 61 Unforested (remnant) Glen Cluanie
## 56 70 Unforested (remnant) Glen Cluanie
## 57 72 Unforested (remnant) Glen Cluanie
## 58 75 Unforested (remnant) Glen Cluanie
## 59 76 Unforested (remnant) Glen Cluanie
## 60 78 Unforested (remnant) Glen Cluanie
## 61 79 Unforested (remnant) Glen Cluanie
## 62 32 Unforested (remnant) Glen Loyne
## 63 33 Unforested (remnant) Glen Loyne
## 64 34 Unforested (remnant) Glen Loyne
## 65 46 Unforested (remnant) Glen Loyne
## 66 49 Unforested (remnant) Glen Loyne
## 67 51 Unforested (remnant) Glen Loyne
## 68 52 Unforested (remnant) Glen Loyne
## 69 1 Unforested (remnant) Glen Mallie
## 70 2 Unforested (remnant) Glen Mallie
## 71 5 Unforested (remnant) Glen Mallie
## 72 15 Unforested (remnant) Glen Mallie
## 73 22 Unforested (remnant) Glen Mallie
## 74 27 Unforested (remnant) Glen Malliepitfall_env$Type <- factor(pitfall_env$Type)
scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_pitfall, type = "n")
with(pitfall_env, points(NMDS_pitfall, display = "sites", col = colvec[Type],
scaling = scl, pch = 21, bg = colvec[Type]))## Warning in plot.xy(xy.coords(x, y), type = type, ...): "scaling" is not a
## graphical parameterwith(pitfall_env, legend("topleft", legend = levels(Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec))
NMDS_pitfall.fit <- envfit(NMDS_pitfall ~ Type , data = pitfall_env, perm = 999)
NMDS_pitfall.fit##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## TypeMature 0.0521 -0.1043
## TypeRegen -0.0412 -0.0321
## TypeUnforested (heath) -0.0783 0.3440
## TypeUnforested (remnant) 0.0355 -0.0632
##
## Goodness of fit:
## r2 Pr(>r)
## Type 0.0528 0.243
## Permutation: free
## Number of permutations: 999plot(NMDS_pitfall.fit, dis = "site")## Warning in text.default(x$factors$centroids[, choices, drop = FALSE], labs$f, :
## "dis" is not a graphical parameterplot(NMDS_pitfall.fit)
### Ground invertebrate NMDS
# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)
pit_ordi_ellipses <- ordiellipse(NMDS_pitfall, pitfall_env$Type, kind = "se", conf = 0.95, draw = "none")
# 2. Extract ellipse coordinates into a data frame
pit_ellipse_df <- do.call(rbind, lapply(names(pit_ordi_ellipses), function(group) {
pit_ellipse_points <- vegan:::veganCovEllipse(pit_ordi_ellipses[[group]]$cov,
pit_ordi_ellipses[[group]]$center,
pit_ordi_ellipses[[group]]$scale)
data.frame(pit_ellipse_points, Type = group)
}))
colnames(pit_ellipse_df)[1:2] <- c("NMDS1", "NMDS2")
library(ggplot2)
# Your main NMDS points
pit_nmds_df <- data.frame(
NMDS1 = NMDS_pitfall$points[, 1],
NMDS2 = NMDS_pitfall$points[, 2],
Type = pitfall_env$Type
)
pit_anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(pitfall_wide)), pitfall_env$Type, distance = "bray", permutations = 9999)
pit_R_val <- round(pit_anosim_plant_com_r_tr$statistic, 3)
pit_p_val <- ifelse(pit_anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(pit_anosim_plant_com_r_tr$signif, 2)))
pit_anosim_label <- paste0("ANOSIM R = ", pit_R_val, "\n", "p ", pit_p_val)
# Site scores
pit_site_scores <- as.data.frame(NMDS_pitfall$points)
colnames(pit_site_scores)[1:2] <- c("NMDS1", "NMDS2")
pit_site_scores$Type <- pitfall_env$Type # Add grouping
# Species scores (usually in $species, or use scores())
pit_species_scores <- as.data.frame(scores(NMDS_pitfall, display = "species"))
pit_species_scores$Species <- rownames(pit_species_scores)
pit_nmds_hull_graph <- ggplot(pit_nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
geom_point(shape = 21, size = 1, stroke = 1) +
geom_polygon(data = pit_ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
alpha = 0.4, linewidth = 0.5) +
geom_text(data = pit_species_scores, aes(x = NMDS1, y = NMDS2, label = Species),
color = "darkcyan", alpha = 0,
inherit.aes = FALSE,
size = 2) +
scale_color_manual(name = "Forest Treatment",
values = colvec) +
scale_fill_manual(name = "Forest Treatment",
values = colvec) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
strip.background = element_rect(fill = "grey", color = "black", size = 1), # Grey background for facet labels
strip.text = element_text(color = "black", size = 12),
panel.spacing = unit(1, "lines"), # Add spacing between panels if needed
panel.background = element_rect(fill = "white", color = "black"), # Apply white background with black border to panels
plot.margin = margin(10, 10, 10, 10) ) +
annotate("text", x = Inf, y = Inf, label = pit_anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
labs(title = "Ground Invertebrate community separation")
pit_nmds_hull_graph + coord_fixed()
save_plot(pit_nmds_hull_graph, "pit_nmds_hull_graph")## png
## 2Ground invertebrate resampling
# Start by calculating the base ANOSIM R-value (with all species included)
base_R_pitfall <- anosim(wisconsin(sqrt(pitfall_wide)), pitfall_env$Type, distance = "bray", permutations = 999)$statistic
# Initialize a data frame to store the jackknife results
jackknife_results_pitfall <- data.frame(Species = character(),
R_value = numeric(),
Delta_R = numeric(),
stringsAsFactors = FALSE)
# Loop over all species (columns in your wide matrix)
species_names_pitfall <- colnames(pitfall_wide)
for (sp in species_names_pitfall) {
# Remove the species (column)
reduced_matrix_pitfall <- pitfall_wide[, setdiff(species_names_pitfall, sp)]
# Remove rows that are now completely zero
nonzero_rows_pitfall <- rowSums(reduced_matrix_pitfall) > 0
reduced_matrix_pitfall <- reduced_matrix_pitfall[nonzero_rows_pitfall, ]
reduced_env_pitfall <- pitfall_env[nonzero_rows_pitfall, ]
# Skip iteration if data becomes too small or empty
if (nrow(reduced_matrix_pitfall) < 2 || anyNA(reduced_matrix_pitfall)) next
# Run ANOSIM on the reduced matrix
anosim_result_pitfall <- anosim(wisconsin(sqrt(reduced_matrix_pitfall)), reduced_env_pitfall$Type, distance = "bray", permutations = 999)
# Store R and change from base
R_val_pitfall <- anosim_result_pitfall$statistic
delta_R_pitfall <- base_R_pitfall - R_val_pitfall
jackknife_results_pitfall <- rbind(jackknife_results_pitfall, data.frame(Species = sp, R_value = R_val_pitfall, Delta_R = delta_R_pitfall))
}
# Sort by greatest impact (largest Delta_R)
jackknife_results_pitfall <- jackknife_results_pitfall[order(-jackknife_results_pitfall$Delta_R), ]
# Plot Jackknife Results
pitfall_jk <- ggplot(jackknife_results_pitfall, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
geom_col(fill = "#0072B2") +
coord_flip() +
labs(title = "Ground invertebrate Order impact on NMDS Clustering",
x = "Species Removed",
y = "Decrease in ANOSIM R") +
theme_minimal()
pitfall_jk
save_plot(pitfall_jk, "pitfall_jk")## png
## 2